Accelerated boundary integral method for multiphase flow in non-periodic 

geometries 



Amit Kumar, Michael D. Graham* 

Department of Chemical and Biological Engineering, University of Wisconsin-Madison, Madison, Wisconsin 53706, USA 



Abstract 

An accelerated boundary integral method for Stokes flow of a suspension of deformable particles is presented 
for an arbitrary domain and implemented for the important case of a planar slit geometry. The computational 
complexity of the algorithm scales as 0{N) or 0(iVlog A^), where A^ is proportional to the product of number 
of particles and the number of elements employed to discretize the particle. This technique is enabled by 
the use of an alternative boundary integral formulation in which the velocity field is expressed in terms of 
a single layer integral alone, even in problems with non-matched viscosities. The density of the single layer 
integral is obtained from a Fredholm integral equation of the second kind involving the double layer integral. 
Acceleration in this implementation is provided by the use of General Geometry Ewald-like method (GGEM) 
for computing the velocity and stress fields driven by a set of point forces in the geometry of interest. For 
the particular case of the slit geometry, a Fourier- Chebyshev spectral discretization of GGEM is developed. 
Efficient implementations employing the GGEM methodology are presented for the resulting single and 
the double layer integrals. The implementation is validated with test problems on the velocity of rigid 
particles and drops between parallel walls in pressure driven fiow, the Taylor deformation parameter of 
capsules in simple shear flow and the particle trajectory in pair collisions of capsules in simple shear flow. 
The computational complexity of the algorithm is verified with results from several large scale multiparticle 
simulations. 

Keywords: accelerated boundary integral, confined, slit, non-periodic, capsule, red blood cells, 
microfiuidics 



1. Introduction 

Multiphase flow in confined geometries is ubiquitous in nature and technological applications. A very 
common example is blood flow in the microcirculation. Recall that blood is primarily a suspension of red 
blood cells (RBCs) in plasma, with the volume fraction (/) of RBCs (hematocrit) typically ranging between 
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(j) ^ 0.1 — 0.3 in the capillaries and reaching as high as i/i w 0.5 in large arteries [l^. The diameter of the 
blood vessels in the microcirculation, which includes the capillaries, arterioles and venules, is typically in 
the range 10 — 125 /iw [l^], such that a discoidal RBC with a typical diameter and thickness of 8 fim and 
2 fim respectively can be strongly to moderately confined. Therefore, any realistic computational study of 
blood flow in the capillaries must account for confinement. Other examples of technological interest where 
confinement effects are usually significant include multiphase fiows in microfluidic devices 5^. Again, 
any realistic model must account for confinement in such problems. Given the importance of multiphase 
flows under confinement, or more generally speaking in non-periodic geometries, it is imperative to develop 
efficient and accurate computational techniques which faithfully represent the system under study, including 
the aspect of system size (meaning number of particles here). The algorithm presented herein has been 
motivated by our goal to study the class of problems described above. We next discuss some related 
previous efforts on the computational studies of multiphase fiows under confinement. 

Boundary integral based methods have emerged as a powerful tool for studying the flow behavior of 
multiphase systems in the limit of negligible Reynolds number, i.e. under Stokes flow conditions. Such 
methods have been employed in the past to study the flow behavior of a variety of particle types including 
drops, capsules, RBCs, and vesicles among others. Most of these prior implementations scale as 0{N^), 
where N is proportional to the number of degrees of freedom in the system. For a system with Np particles, 
each of which have been discrctized into A^a elements, the number of degrees of freedom in the system 
scales as TV ^ NpN^- The 0{N^) scaling above assumes an iterative solution of the discrctized system of 
equations, where the number of iterations is independent of iV; a direct solution will result in a scaling of 
0{N^), while a system size dependent number of iterations with an iterative solution results in a scaling 
higher than O(A^); the worst case scaling being 0{N^). The 0{N'^) scaling is usually prohibitive, such that 
it precludes a numerical study of large system sizes. It is therefore not surprising that many of the past 
studies have been limited to an 0(1) number of particles. 

To overcome these limitations, there have been several efforts to develop accelerated techniques, where an 
accelerated techniciuc is assumed to give a scaling closer to the ideal 0{N), while being sufficiently accurate 
at the same time. These accelerated techniques employ either some variant of the particle-particle-particle- 
mesh (P^M) method ^, or the Fast Muflipolc Method (FMM) One of the earliest Stokes flow boundary 
integral implementation with acceleration was perhaps presented by Greengard et al. [l8|, who employed 
the FMM for acceleration in complex domains. Using the particle- mesh-Ewald (PME) method, Metsi [46 



developed an accelerated implementation of the Stokes flo w boundary integral m ethod for her studies on two 



dimensional periodic suspensions of emulsions and foams. 



Zinchenko and Davis 



6J, |65[ employed multipole 



expansion accelerated boundary integral method to study large number of drops in a periodic geometry 
under shear. Freund {igI ] used the smooth particle-mcsh-Ewald method to study the motion of periodic 
suspensions of RBCs and leukocytes in two dimensions; this was later extended to three dimensions by Zhao 

2 



et al. 



63|. In the latter study 



63|, the effect of confinement was incorporated by explicitly discretizing 
the walls, which generally has unknown tractions and known no-slip velocity conditions. This explicit 
discretization is necessary because the periodic Green's function docs not inherently satisfy the no-slip 
condition on the walls. Additionally, the previous authors employed a staggered time integrator, such that 
the wall tractions and the particle surface velocities were not determined simultaneously; this is due to 
the large cost associated with their simultaneous solution. Another potential drawback with the periodic 
Green's function is that it has a zero mean flow and a non-zero mean pressure gradient associated with it 



22l |. As a consequence, the pressure drop in the system is not directly a specified quantity, and must be 
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Note that many experiments 



solved for by varying the mean flow, which is a specified quantity 
on pressure driven flow have a specified pressure drop, and it is therefore desired to specify the pressure 
drop directly in numerical simulations without incurring additional computational costs. We also remark 



that the specified mean flow includes the fiow outside the walls, as that is a part of the simulation box |63j : 
consequently, in such a method, neither the mean flow between the walls nor the pressure drop is a directly 
controlled quantity. Other recent work of possible interest is Rahimian et al. [s^ , where a FMM accelerated 
boundary integral method is presented. While this implementation was developed for an arbitrary domain, 
its applicability is restricted to two dimensional systems. In a subsequent article |62l |. the previous authors 



generalized their implementation to three dimensions, though only an unbounded domain was considered. 
At this point, it must be emphasized that all the prior accelerated implementations of the boundary integral 
method are based on either the free space Green's function or the periodic Green's function; in such a case, 
the boundaries of the confined domain arc required to be explicitly discretized. 

We next discuss previous boundary integral implementations employing the Green's function for the 
geometry of interest. Such a Green's function satisfies the appropriate boundary conditions at the domain 
boundaries; consequently, the unknowns at the domain boundaries, e.g. hydrodynamic tractions, do not 
enter the boundary integral equation. A popular geometry for which several boundary integral implemen- 
tations have been developed is a slit - the region between two parallel walls. The Green's function for this 
geometry has been provided by Liron and Mochon 42| . A boundary integral implementation based on this 



58| for rigid particles. This was later extended b y G riggs 



Green's function was developed by Staben et al 
et al. 21[ for studies on drops in the same geometry. In a related work, Janssen and Anderson 



271 also 



implemented a boundary integral method for drops between two parallel walls, though that was restricted to 
matched continuous and dispersed phase viscosities. This was later extended by Janssen and Anderson [28 
to include non-matched viscosity problems. It is important to emphasize that none of these implementations 
are accelerated and have a computational cost of at least 0{N'^). Consequently, it is not surprising that 
all of the studies described above were limited to a few particles and are thus not suitable for studying 
suspension dynamics. 

We now briefiy discuss examples of other numerical techniques employed in the literature for studies on 
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the flow behavior of particles under confinement. One such simulation technique is the immersed boundary 
method, which has been employed, e.g. by Doddi and Bagchi [ll| and Pranay et al. jsij, for studies on 
capsules in a slit. Another popular technique is the lattice-Boltzmann method. As an example, MacMeccan 
et al. [4J] developed a coupled lattice-Boltzmann and finite element method to study deformable particles, 
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which included studies under confinement. A somewhat related algorithm is the multiparticle collision 

to study RBCs and vesicles in capillaries. 



dynamics, which has been used by, e.g., Noguchi and Gompper 

The ideal computational cost of all the above numerical techniques is 0{N). Finally, Swan and Brady 
have developed an accelerated Stokesian dynamics method for rigid spherical particles in a slit that uses 
ideas related to those presented here and scales as 0(A^log A^). 

In the present study, we develop an accelerated boundary integral method for multiphase flow in an 
arbitrary geometry and implement it for a slit geometry as shown in Fig. [TJ The computational complexity 
of this algorithm scales as 0{N) or O(A^logA^) depending on the specific numerical scheme employed. 
The latter scaling of O(A^logA^) is associated with the use of fast Fourier transforms (FFTs) if one or 
more directions have periodic boundary conditions, though that is not a requirement of our method. In the 
present effort, we provide a detailed description using the example of the slit geometry; its extension to other 
geometries is straightforward. The acceleration in our method is provided by the use of General Geometry 
Ewald-Like Method (GGEM) [2^. The choice of GGEM as the acceleration technique necessitates the use 
of an alternative boundary integral formulation in which the velocity field is expressed solely in terms of a 
single layer integral jsol l ; its unknown density is obtained from a second kind integral equation involving the 
double layer integral. The resulting single and double layer integrals are computed efficiently employing the 
GGEM methodology, which results in the aforementioned favorable scaling of 0{N) or 0(A^log A^). 

The organization of this article is as follows. In Sec. ([5]), we provide a brief overview of the GGEM ac- 
celerated boundary integral method and also discuss some of the limitations of the current implementation. 
Following this, in Sec. we present the boundary integral formulation and discuss its numerical imple- 
mentation using GGEM. In Sec. we present the procedure to compute the hydrodynamic traction jump 
at a particle surface using the example of a capsule with a nco-Hookean membrane. The solution procedure 
for the discretized boundary integral equation is presented next in Sec. (O. An extensive validation of our 
method is presented in Sec. (jO]). Lastly, in Sec. we present results from several large scale multiparticle 
simulations and verify the computational complexity of our algorithm. 



2. Overview of the current work 

We summarize here some of the key aspects of the present work. Each of these points are discussed 
further later in the article. 

• The general geometry Ewald like method (GGEM) employed in the current work for acceleration 
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essentially yields the geometry- dependent Green's function and other associated quantities. This work 
is a first instance of an accelerated boundary integral method based on the geometry-dependent Green's 
function. Prior implementations have employed either the free-space Green's function (in case of FMM 
accelerated methods) or the periodic Green's function (in case of PME accelerated methods). 

The GGEM methodology decomposes the overall problem into a local problem and a global problem, 
essentially by splitting the Green's function into local (singular but exponentially-decaying) and global 
(smooth but long-ranged) parts. The implementation of the local problem is similar to that of the 
traditional boundary integral method. However, since the local Green's function decays exponentially 
with distance from the source of the singularity, distant elements are not coupled and the local solution 
can be obtained in 0{N) operations. 

The global problem involves solving a single phase Stokes equation in the domain of interest with 
known boundary conditions and with a known smooth distribution of force densities. It is in this 
problem that the coupling between distant elements appears. In solving the global problem, one is not 
concerned with the particle interfaces, or the different viscosity fluids present inside and outside the 
particle (if that is the case in the original problem) . This major simplification allows the use of a wide 
variety of fast and accurate numerical techniques present in the literature for the solution of Navier- 
Stokes equation in an arbitrary domain. All these methods, including those based on finite difference, 
finite volume, finite element, and spectral methods are suitable here [6, tI-Ii^Ii^. In addition, various 
fast and efficient implementations of Navier-Stokes solvers on GPUs and distributed memory systems 
are readily available 
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• We present a spectral 0(A^log A^) Stokes flow solver for the global problem in a slit geometry, which, 
as indicated in the introduction, is one of the most widely-studied confined geometries studied in 
the literature. This solver employs a Fourier-Chebyshev Galerkin method in conjunction with the 
influence matrix approach [3]. The unknown coefficients of the Fourier-Chebyshev series expansion 
are computed with a direct 0{N) algorithm - no iterations are necessary here, as would be the case 
with FMM or PME accelerated methods in a slit. 

The implementation of the methodology presented in this paper also has some limitations. An important 
limitation of the current implementation concerns the evaluation of the near singular integrals - these 
integrals arise when the gaps between the particles become very small, and, if not treated appropriately, 
may cause the simulations to diverge. In the present work this is alleviated by requiring that the minimum 
interparticle gap in the system be always maintained above a specified value; this is achieved by the use of 
an overlap correction procedure in an auxiliary step. There are several other minor limitations of the current 
implementation. For example, we currently use linear elements to discretize the particle surface. It may be 
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Figure 1: Schematic of the problems considered here: a dispersed phase with viscosity /i inside the domain boundaries denoted 
by , containing (for example) two particles with internal viscosities Ai/i and A2/i respectively; their surfaces are denoted by 
Si and 52- The undisturbed flow is denoted by u°°. 



beneficial to employ higher order discretizations, like a spectral discretization, which could be particularly 
helpful for the accurate evaluation of the near singular integrals. In this paper results are reported only for 
a slit geometry. It will be appropriate to develop efficient implementations of the our methodology for other 
geometries like a cylinder. It must be emphasized that none of the above limitations are inherent to our 
methodology and we hope to address these in future efforts. 

3. Problem Formulation and Implementation 

3.1. Boundary integral equation for fluid motion 

We consider a three-dimensional suspension of deformable particles (e.g. fluid-filled capsules) as shown 
in Fig. ([T]), where both the suspending fluid and the fluid enclosed by the particles are assumed to be 
Newtonian and incompressible. The viscosity of the suspending fluid is taken to be fi, while the viscosity of 
fluid enclosed by capsule m is taken to be Am/i, such that A,„ is the viscosity ratio of the interior and the 
exterior fluid for this particular capsule. The Reynolds number for the problem is assumed to be sufficiently 
small that the fluid motion is governed by the Stokes equation. Under these assumptions, one may write the 
velocity at any point in the domain with an integral expression involving only the boundary of the particles 
[soi l. We flrst introduce the formulation that is most commonly used, 

2 1 
«,(xo) = — — (xo) - — — — — ^ / A/,(x)G,,(xo,x)d5(x) 

(1) 

+ TTTTItI^^^"^") / "i(x)Tyfc(x,xo)nfc(x)d5(x), 
47r(l -I- A„,) Js^ 

where u(xo) is the fluid velocity at a point xq lying on the boundary of particle m (i.e. xg E S"^, S"^ 
denotes the surface of particle m). u°°(xo) is the undisturbed fluid velocity at the point xq, Af(x) is the 
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hydrodynamic traction jump across the interface jsot . and the sums arc over all the Np particles in the 
system. The Green's function and its associated stress tensor are denoted by G and T respectively in the 
above equation, and integrals involving them as the kernel are typically referred to as the single layer integral 
and the double layer integral respectively jsol. Isojl. From here onwards, a principal value of the double layer 
integral over a part of the boundary is assumed whenever the target point xq lies on that boundary. For 
example, in the above equation, the double layer integral over S"" is assumed to denote the principal value 
when n = TO. A crucial aspect of the above formulation is that the Green's function G and associated stress 
tensor T are taken to satisfy the boundary conditions imposed at the system boundaries, so the integrals 
above only involve the internal (intcrfacial) boundaries; if the Green's function for any other geometry is 
employed (e.g. periodic), additional integrals over the domain boundaries arise in Eq. ([T]). 

The above form of the boundary integral equation ([1]), using the free space Green's function (the Oseen- 
Burgers tensor) or the Green's function for a triply periodic domain given by Hasimoto 2^ is widely used 
in the literature and is the basis for numerous numerical implementations, including the references cited in 
the introduction. However, for reasons that will be discussed shortly, this form is not amenable to numerical 
solution by an accelerated method in an arbitrary domain when using the Green's function for that domain. 
In the present effort, therefore, we employ an alternative formulation in which the fluid velocity is expressed 
solely in terms of the single layer integral with density q(x) as follows: 



Uj(Xo) = 



(xo) + V / g.(x)G,,(xo,x)d5(x). (2) 



The single layer density q(xo) satisfies (for xq G 5*'" ) 



<lj i^o) + X^"-fc(xo)X! / '7j(x)Tj^fe(xo,x)d5(x) = -- — 



A/,(xo) 
A™ + 1 



where k,„ is defined as 



K„J-(xo) , (3) 



(4) 



while f°° is the traction at a given point (computed with the suspending fiuid viscosity n) due to the 
stress generated in the fiuid corresponding to the undisturbed fiow u°° (see [Appendix B] for examples). 
In lAppendix A| this formulation is derived from the Lorentz reciprocal theorem for the case of a single 

. . . n . . 

particle. This derivation follows closely the approach outlined in Chap. (5) of Pozrikidis |50j and is provided 
here for completeness. 

We now clarify the motivation for employing Eqs. ^ and ^ rather than the more commonly employed 
formulation in Eq. ([T]). We begin by noting that the first argument xi of G(xi,X2) and T(xi,X2) denotes 
the field (target) point of the functions, while the second argument X2 denotes the location of the pole 
(source) of the singularity that drives the fiow. A close look at Eq. ([T]) reveals that the operand of G(xo, x). 
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A/i(x), is a function of the position of the pole of the singularity (x) and that the field point of the G tensor 
is same as the target point of the overall boundary integral equation (xq). In other words, the operand of 
G is independent of its target point, and consequently the same collection of point forces can be used to 
compute the velocity at any target point. This requirement is essential to any accelerated method, as in 
such methods a part of the calculation gives the velocity (or other relevant quantities) simultaneously at 
all target points (e.g., boundary element nodes) due to all the singularities present in the system. This is 
possible only if the operands of the singularities are independent of the target points of the singularities, 
and instead are functions of the location of the pole of the respective singularities. 

With this aspect clarified, it is seen that this important condition is not satisfied for the double layer 
kernel T(x, xg) in Eq. ([ij as its multiplicands u(x) and n(x) are functions of its target point x. Also note 
that no general relationship exists that would allow one to switch the location of the pole and the field 
points in T. (This is possible for G, since, by self-adjointness, Gij(x, xq) = Gji(xo,x) {s^])- Hence the 
above formulation ([1]) is not suitable for our purposes here, though it can still be used for problems in which 
the viscosity ratio is unity, as the double layer integral vanishes in such a case [si! ]. In contrast, in the 
formulation employed in this work (Eqs. [2]and|31), the multiplicand q(x) of both G(xo,x) and T(xo,x) is a 
function of the location of the source point x. Hence, it is amenable to numerical solution by an accelerated 
method. 

We now describe the fast computation of the velocity and pressure fields due to a collection of known 
point forces, which is closely related to the problem of computing the Green's function and its associated 
stress tensor in the geometry of interest. Later in Sees. (j3.4p and p.5p . we employ this technique to compute 
the single layer and double layer integrals. 

3.2. GGEM Stokes flow solver for a collection of point forces 

Consider the velocity field u(x) and the pressure field p(x) due to a collection of Ng point forces, such 
that the strength and location of the v*^ point force are given by and x" respectively. The velocity and 
pressure fields above are obtained from the solution of the Stokes and the continuity equation as shown 
below: 



and subject to given boundary conditions on the system boundary S . By definition, the above velocity 
and pressure fields along with the associated stress tensor cr can be written in terms of a Green's function 




(5a) 



V • u(x) = 0, 



(5b) 
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Figure 2: Variation of the global Pg(x) and the local pi(x) force density along the x-axis, given the center of the force density 
is at the origin. Note that both of these densities are functions only of the distance r from the origin (see Eq. [S]!. Also note 
that Pg(r) + Pi(r) = i5(r). In numerical calculations, we set pg{r) = for r > A/a. For plotting Pg{x) and pi{x) here, we set 
a = 1. 



G, its pressure vector P and stress tensor T as 



N 



p(x) 



(6a) 
(6b) 



i/=i 



cr,fc(x) ^ ^ X]^y'^(^'^'')^J■ 



The stress tensor Tijk in the above equation is obtained from Gij and Pj from the Newtonian constitutive 
equation 

T,-,(x,x ) = -P,(x,x g,^, + j- (7) 

A close look at the boundary integral equations ^ and ^ shows that to evaluate the integrals we do not 
explicitly need the Green's function G and its stress tensor T but only their products with the density q. 
Put simply, our end goal is to quickly find in time 0{Ns) or 0{Ns\ogNs) the velocity u and the stress 
tensor cr due to a given set of point forces - explicit construction of G and T are not necessary. 

One of the attractive features of the method presented here lies in the fact that it is applicable to an 
arbitrary geometry. For simplicity and considering the interest of the present work, we provide a detailed 
discussion only for a slit geometry (see Fig. [S]); generalization of the formalism for an arbitrary geometry 
is straightforward 



25| . For the present slit domain, there is a no slip boundary condition at the two rigid 
walls at y = and y ~ H, while periodic boundary conditions are assumed in the other two directions x 
and z, with spatial periods and L^, respectively. 

To achieve the computational complexity of 0{Ns log A^^) alluded to above, we employ the general geom- 
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etry Ewald like method (GGEM) f24j for computing the velocity and stress fields due to a given collection 
of point forces. We briefly describe GGEM next. In the GGEM methodology, the Dirae-delta density in 
Eq. ([5]) is expressed as the sum of a smoothly varying quasi- Gaussian global density Pg(f) characterized by 
a "splitting parameter" a and a second local density pi{r) (see Fig. Here f is a position vector relative 
to the pole of the singularity, f = x — x"^ . The above global and local densities are respectively given by the 
following expressions: 



pKf)-<^(f)-Ps(f), (8b) 

where represents a length scale over which the delta- function density has been smeared using the quasi- 
Gaussian form above, and consequently it also represents the length scale beyond which both the global 
and the local densities are effectively zero. It is important to emphasize that the total density remains a 
(5- function, i.e. Pg{r) + Pi{y) = 5{r). The motivation for this particular splitting of the (5-function density 
into Pg{v) and p;(f) will be obvious below. 

We next consider the solution of the Stokes and continuity equation with the above two force densities 
as forcing functions. The solution driven by the local density, u'(x). p'(x), and (t'(x) (velocity, pressure, 
and stress respectively) will be referred to as the local solution, and satisfies the local problem 

- Vp'(x) + ^V2u'(x) + ^ gVz(x - x"^) = 0, (9a) 

V-u'(x) = 0. (9b) 

This equation will be solved in an unbounded domain, i.e. the solution decays to zero at infinity. The 
solution u^(x), p^{x), and (x) driven by the global density will be referred to as the global solution, and 
satisfies the global problem 

- Vp9(x)+^V2u9(x) + 2gVg(x-x'^) =0, (10a) 

i/=i 

V-uS(x) = 0. (10b) 

The boundary conditions for the global problem are set so that the total velocity field u(x) = u'(x) + (x) 
satisfies the specified boundary conditions for the overall problem. Once the local and the global solutions 
are known, the solution to the overall problem is obtained as 

u(x) = u'(x) + uff(x), (11a) 

p(x)=y(x)+p9(x), (lib) 

(t(x) = (t'(x) +cr9(x). (11c) 
10 



We next discuss the solution procedures for the local and the global problems. 
3.2.1. Local solution 

Consider first the local problem. The solution to this problem, u'(x), p'(x), and c'(x) is expressed by 
a set of equations similar to that in Eqs ([5]), which, for the simplicity of nomenclature, is called the local 
Green's function G' and its associated quantities. In short, we append the superscript / to the previously 
defined quantities to denote the solution associated with the local density as the forcing function, and these 
are given by the following: 



GL(x,x^) = + ^) erfc(af) - ^ [s^, - ^) 



, 2 -.2 



e"" , (12a) 



P'(x,x^) = ^erfc(af) + 1^ ( 1 - ) e""'^', (12b) 



ml / u\ QXiXjXji , ^, 12q; XiXjXf^ _ 2-2 

4Jx,x-) = ^erfc(ar) - -^^^e 



(12c) 



Z' c - , c- J e 2,XiXjXk I 

— ^ I OjkXi + dikXj + dijXk ) e 



2 "2 

■a r 



where x = x — x'^, while r ~ |x|. The velocity and stress fields are then obtained as: 



v=l 



aU^)^^Y^T^,,{^,^n9^. (13b) 

The solution in Eq. (|12p has been obtained with free-space boundary conditions, i.e. all of them decay to 
zero at infinity. In other words, the local solution is independent of the geometry of interest. The violation 
of the boundary condition requirements of the domain by employing free space boundary conditions above 
will be corrected by appropriately choosing the boundary conditions for the flow problem associated with 
the global force densities as the forcing function. 

An important observation at this point is that the local solutions in (|12p are short ranged, decaying 
approximately as e^" . Consequently, the contribution from the local solution can be neglected beyond a 
length scale ^ from the origin of the corresponding local density. In this work, this cutoff length was 
taken as Vcut = 4/a throughout. The near neighbor list required for the efficient computation of the local 



solution is generated by the 0{Ns) cell-linked list algorithm 

It is important to point out that the G' in Eq. (|12b,) has the same functional form as the real space 



term in the periodic Stokeslet (Green's function) provided by Hasimoto [22|. In other words, Hasimoto's 



solution for the periodic Stokeslet can also be obtained by first splitting the (5-function density into the local 
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and global densities as in Eq. ([S]); the local problem is then solved as described above, while the global 
jroblem is solved with a Fourier Galerkin method with the appropriate assumptions described in Hasimoto 



221 ■ Since PME accelerated methods (e.g. [3j|) for Stokes flow employ the periodic Stokeslet given by 



Hasimoto [63|, this observation illustrates a connection between PME like methods and GGEM. A very 
important distinction, though, is that the performance of PME like methods is tied to the use of discrete 
Fourier transforms and thus periodic domains, which is not the case with GGEM (discussed below) and 
hence the latter's much broader applicability. 

3.2.2. Global solution 

We now describe the solution to the global problem, i.e. the flow problem associated with the collection 
of global force densities. We first discuss the boundary conditions for the global problem. As was mentioned 
earlier, the overall solution for a given collection of point forces is the sum of the corresponding quantities 
from the local and the global solutions, see Eq. ([TT]) . It is obvious that the same should be true for boundary 
conditions. Consequently, to satisfy any type of boundary condition (e.g. Dirichlet) at an arbitrary location, 
we set the boundary condition for the global part so that its sum with the known contribution from the 
local part (above) adds up to the required value. Again, we employ the example of the slit geometry, noting 
that this scheme is equally applicable to other geometries. To satisfy the no-slip condition at the two rigid 
walls of the slit, we require the following at y = and y ~ H: 

uf = -u'. (14) 

Note that in the present formulation the static no-slip condition is always imposed at the rigid walls for 
computing the velocity field due to the Green's function (or point forces). This is true even in problems where 
the walls may not be at rest, a common example being simple shear flow. The effect of the undisturbed 
flow enters the boundary integral equation via u°° and f°° in Eqs. ^ and ([3]). To satisfy the periodic 
boundary conditions in x and z directions, we impose equivalent periodic boundary conditions in the global 
calculation. As far as the local solution is concerned, we require that it decays to a negligible value over 
a length scale equal to half of the spatial period in x and z directions or smaller: i.e. Vcut < and 
Tcut < Lz/2. Given the above choice rcut = 4/a, we require that aL^ > 8, aL^ > 8. This fact, coupled with 
the minimum image convention Q employed in the computation of the local solution ensures its periodicity. 

Before proceeding further, we note that, in general, two subtleties arise in considering the behavior of 
the global solution near boundaries. The first is the issue of boundary shape. In the present work we take 
the boundary to be smooth on the scale of the suspended particle size. If that is not the case, it will be 
necessary to resolve the length scales of the boundary roughness. For such a boundary it might be convenient 
to revert to a conventional accelerated method and explicitly discretize the boundary. Alternately, in the 
present context, the global solution could be obtained using a locally-refined mesh near the domain boundary 
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to capture its features (see, e.g. Fard et al. without destroying the seahng of the method with the 

number of suspended particles. 

Another issue arises in principle if a particle very closely approaches a (smooth) boundary. This does not 
arise in the context of the present application as deformable particles migrate away from solid surfaces in 
shear flow due to the hydrodynamic dipole interaction between the particle and the wall 157|. Nevertheless, 



it can happen in principle and leads to the situation where the boundary condition for the global problem 
that must be satisfied becomes nearly singular. This is because, at a no-slip boundary, we require = — u' 
and a point force a distance e away from the boundary leads to a local velocity u' of 0(l/e) on the boundary. 
This situation can be addressed by adding the image system for a plane wall [sj (and splitting it into local 
and global parts), in which case the effects of the singularities cancel on the wall. See Swan and Brady 
for a related discussion in the Stokesian dynamics context. 

Having discussed the boundary conditions for the global problem, we turn to the solution procedure 
of the Stokes and the continuity equation with the given collection of global force densities as the forcing 
function; see Eq. (|10|) . For an arbitrary geometry one may employ any desired discretization scheme 
for the solution of the global problem. If a finite difference or a finite element scheme is used, then the 
solution can be obtained at a cost of 0{N) when the resulting sparse matrix equations are solved iterati yely 
with proper preconditioners; the multigrid preconditioner for Stokes flow is an attractive choice [l^ IsgI ]. 
Section 17.21 contains further discussion of the scaling of computation time with problem size. For the slit 

employed discrete Fourier series approximation in the periodic 



problem of interest here, past work 
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X and z directions, while a second order finite difference discretization was employed in the wall normal y 
direction. In the present work, we develop a fully spectral solution procedure by employing the discrete 
Chebyshev polynomial approximation ^ in the wall normal direction, while the discrete Fourier series 
approximation is used in the periodic x and z directions. For example, the x-component of the global 
velocity = {u^ ,w^) in Eq. ([TU)) is expressed as 

Ny-1 

"'W= E E E i^LnTniy)e^'^'^^^^e^'''"'^/^^, (15) 

l = -N^/2 m=~N^/2 n=0 



where T„(j7) = cos(ncos~^ y) is the Chebyshev polynomial of the n*"^ degree [6|, y represents the mapping 
from [0,i7] to [—1,1]: y = 2y / H — 1, while N^., Ny, and respectively denote the number of terms 
(modes) in the corresponding series approximation. Similar expressions are written for other components of 
the velocity and the pressure. An important implication of this representation, particularly with regard to 
the pressure, is that the pressure drop associated with this point force solution is always zero over the spatial 
period of the domain, while the mean fiow is (in general) non-zero. This ensures that the pressure drop 
obtained from the boundary integral method always equals the pressure drop specified in the imposed bulk 
flow (i.e. in the absence of the particles). Returning to the expression in Eq. p3|) . we note that the use of 
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the Fourier series approximation in x and z directions ensures that the periodic boundary conditions in these 
directions are inherently satisfied. The Chebyshev polynomials, on the other hand, do not automatically 
satisfy the boundary conditions in the wall normal direction; the satisfaction of these boundary conditions 
was accomplished by employing the tau method In the tau method, the equations for the highest two 



modes in the series approximation are replaced by equations representing the two boundary conditions; see, 
e.g., Canuto et al. Q or Peyret for details. An attractive feature of the discrete Chebyshev polynomial 
approximation is that FFTs can be used for rapidly transferring information from the physical to spectral 
space and vice- versa Q. A major drawback, though, with solving differential equations with Chebyshev 
polynomial approximation is that the differentiation matrix is full in both the spectral and the physical 
space [4^ (in contrast, the Fourier differentiation matrix is diagonal in the spectral space). Due to the 
full nature of the Chebyshev differentiation matrix, a straightforward implementation for solving the Stokes 
flow problem will lead to an 0{Ny) method. For the incompressible Stokes flow problem here, though, an 
alternate approach exists in which the solution to the Stokes equation is obtained from the solution of a series 
of Helmholtz equations In this case, with a little manipulation, a quasi-tridiagonal system of equations 
results, for which a direct 0{Ny) algorithm exists [4^. This approach for solving the incompressible Stokes 
equation, or more generally the incompressible Navier-Stokes equations, is popularly known in the literature 
as the Kleiser-Schumann influence matrix method Q. A detailed discussion of this approach including the 



equations being solved and their respective boundary conditions is presented in Appendix C Here we only 
sketch out the main computational aspects of this approach. To begin, each of variables appearing in the 
Stokes and the continuity equation are expanded in a truncated Fourier series in x and z directions; see, e.g., 
Eq. (|15[). These expressions are then substituted in the Stokes and the continuity equations. Subsequently, 
by the application of the Galerkin method, a set of coupled ordinary differential equations (ODE) in y is 
obtained for each of the Fourier modes of all the unknown variables (velocity components and pressure). 
These coupled ODEs are solved with the Chebyshev-tau influence matrix method, which involve Chebyshev 
transformations, quasi-tridiagonal matrix equation solves, and inverse Chebyshev transformations in that 
order. The solution thus obtained yields the Fourier coefficients of the velocity components and the pressure. 
An inverse Fourier transform then leads to the solution for the velocity and the pressure in the physical space. 
The computation of the stress tensor requires the derivatives of the velocities; these differentiations are 
performed in the transform space [4^ . All of the above Fourier and the Chebyshev transforms along 
with their inverse transforms are performed using the EFT algorithm. Thus, the asymptotic computational 
cost of the solution procedure for the global problem scales as NlogN, where N — N^NyN^. Assuming 
N ^ Ng (see Sec. 17. 2|) . we obtain the asymptotic computational cost of the global solution as 0{Ns log A^s)- 
A further discussion on the computational complexity of the algorithm is provided in Sec. (j7.2p . 

We next introduce some of the important parameters associated with this solution procedure. Associated 
with each of the and Nz Fourier modes, there are and Nz equispaeed trapezoidal quadrature points; 
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(a) Observation point is a mesh point (b) Observation point is not a mesh point 

Figure 3: Relative error in the velocity (Au/u) due to a point force at a given observation point in a slit geometry with = 
H = Lz = L. The abscissa in the plots is (aAj/m)~^, while different curves are for different values of Tent as labeled in the key. 
The strength of the point force is given by (_F, F, F), while its coordinates arc (0.25L, 0.25L, 0.25L). (a) The observation point 
is (0.5-L, 0.5L, 0.5-L), which is always maintained on a mesh point, and (b) the observation point is (0.4312L, 0.3734L, 0.5234L) 
which is not a mesh point. In (a) we also plot the function y ~ e"^"^ (chosen to closely match other curves on the plot) 
to demonstrate the exponential convergence of the solution, while in (b) we plot the function y ~ to demonstrate the 
algebraic convergence due to dominant interpolation errors at higher values of {aAym)~^ ■ 

the corresponding spacings arc denoted by Axm = L^-jN^ and Az,„ — L^/Nz- Similarly, associated with 
the Ny Chcbyshcv polynomials, there arc Ny Chcbyshev Gauss-Lobatto quadrature points, the j*'' of which 
is given by yj = H/2{1 + cos(7r(j — l)/{Ny — 1))); the mean mesh spacing in this case is denoted by 
= H/{Ny — 1). Unless otherwise mentioned, the mean mesh spacings in all three directions are kept 
equal in simulations, i.e. Axm = Ay™ — Az^. For computing any of the above transforms, the value of 
the corresponding physical variable is required only at the corresponding quadrature points. Likewise, as is 
customary, the final solution for the velocity, pressure and stress are also computed only at these quadrature 
points. This last step is essential in maintaining the optimal computational complexity of O(A^logA^) 
alluded to above. The velocity and stress at any point not on the mesh is obtained via interpolation; here 
we employ 4*'' order Lagrange interpolation for which the error decays as ft,^, where h is the characteristic 
mesh spacing. The error, therefore, is expected to decay exponentially fast with the number of modes for 
any point on the mesh, while it is expected to decay as for any point not on the mesh. It is appropriate 
to pointout here that exponential convergence of the solution is possible even at a non-mesh point while 
maintaining the computational complexity of O(A^logiV) - this can be achieved by employing the basic 
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for details. 



principles of non- uniform FFT calculations; see, e.g., 

3.2.3. Convergence of the G GEM solution 

We demonstrate next the convergence behavior of the GGEM Stokes flow solution presented above. It 
will be helpful to begin this section with a discussion of various sources of error in the solution procedure. 
It should be obvious that the overall error in the solution results from errors in both the local and the 
global solution procedures. The error in the local solution arises due to its truncation beyond a distance 
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of Tcut from the source of the singularity - this error scales as e""^*"^"* (Sec. 13.2.11 40|). Typically, we set 
Tcut = 4/q!, which is expected to result in an error of 0(10^^). Smaller error in the local solution can be 
obtained by increasing the value of rcut- The error in the global solution has three different sources. The 
first source of error is the truncation of the Fourier-Chebyshev series expansion at some finite number of 
modes; see Eq. (|15|) . The error due to this truncation is expected to decay exponentially fast on the mesh 
points with the parameter {aAym)~^ - this convergence rate results due to the spectral nature of the global 
solution procedure employed here. For other solution procedures, such as a finite difference scheme |25l. l51j. 
an algebraic convergence with the parameter (aAy„i)^^ will be obtained. Next, for any point not on the 
mesh, the global solution has to be obtained by interpolation from nearby mesh points, which introduces an 
additional error scaling as (aAy„i)^ in the current work. Lastly, there is also an error associated with the 
assignment of the global force density on the mesh points - this is required for computing its Fourier and 
Chebyshev transforms. The global density decays approximately as e"^"''-' with distance r from the origin 
of the density (Eq. Therefore, just like the local solution, we truncate the global density beyond a 

distance of Vcut ~ 4:/a from the origin of the singularity. It is expected that the error due to this truncation 
will scale as e"'"'' . Again, a smaller truncation error can be obtained by setting a larger value of rcut ~ 
in such cases the cost associated with the global force density assignment on the mesh points can be large 



and the fast Gaussian gridding algorithm 
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is recommended. 

Having discussed the various sources of errors in the solution procedure, we turn to verifying the expected 
convergence behavior with a test problem. In this test problem, we compute the velocity field due to 
a point force in a slit of side L, i.e. ~ = H = L. The point force is located at coordinates 
(0.25L, 0.25i, 0.25L), while its strength is given by {F,F,F). The velocity due to this point force will be 
computed at two target points. The first target point is located at the center of the box (0.5L, 0.5i, 0.5L), 
which is easy to maintain on a mesh point, while the second target point is a randomly chosen point with 
coordinates (0.4312L, 0.3734i, 0.5234L), which is unlikely to be a mesh point. In this test study, we will 
keep the value of a fixed at a = 80/3i, while the value of rcut ~ C/a (C is a constant) and the mean mesh 
spacing a Ay™ will be varied. Figure shows the relative error in the x component of the velocity Au/u 
as a function of {aAy„i)^^ for several different values of rcut for the first target point, while the same is 
shown for the second target point in Fig. ([SJd) . The data points in these plots were obtained by varying Ny 
between 17 and 181, while the solution computed with Ny = 225 and rcut ~ 0.5L (i.e., C = 40/3) is taken 
as the reference for computing the relative error. Focusing first on Fig. ([3^), we observe an exponential 
convergence in the velocity with increasing (aAym)^^ , though it eventually levels off at a value depending 
on the choice of rcut- For the typical value of rcut = 4/a, an error of O(10~^) is obtained. We next focus 
on the velocity convergence at the second target point (Fig. Eb), which is a non-mesh point. For the 
choice rcut = 4/a, we again observe an exponential convergence initially with increasing {aAyni)~^ that 
eventually levels off at approximately the same value as for the first target point. For higher values of rcut 
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(a) Triaiigulation 



(b) Parent Triangle 



Figure 4: (a) Discretization of a sphere into triangular elements. The number of triangular elements and the number of vertices 
for the discretization shown in the figure are A^a = 1280 and Ni, = 642 respectively, (b) A schematic of the parent triangle. 
Edges 1-2 and 1-3 arc of unit length. 

(C = 5 and C = 6), we observe an exponential convergence initially, though a convergence rate scaling as 
(aAym)^ is observed at higher values of {aAym)~^ ■ The latter convergence rate results from interpolation 
errors becoming dominant at higher values of (aAy,„)~^. We also note that an exponential convergence 
is observed only when the length scale of the global force density is well resolved by the numerical mesh. 
Since the length scale of the global force density is represented by a^^, the requirement for an exponential 
convergence is quantitatively expressed by the condition {aAy„i)^^ > 1 (i.e. aAy,„ < 1). This requirement 
on the mesh spacing is more easily appreciated if one interprets (aAy„i)~^ as the number of mesh points 
per unit smearing length represented by a~^; therefore the larger is (aAym)^^ , the higher is the resolution 
of the numerical scheme. Based on extensive numerical tests presented in this paper, aAym = 0.5 is a 
recommended value as convergence was usually observed at this resolution. 

3.3. Surface discretization 

Having described the procedure for the fast computation of the velocity and the stress fields associated 
with a given collection of point forces, we now turn to the numerical solution of the boundary integral 
equation introduced in Sec. In this section, wc describe the discretization of the particle's surface 

into elements along with the basis functions employed over each element. Following this, we describe the 
numerical implementation of the single and the double layer integrals present in the boundary integral 
equation. It should be emphasized that accelerated approach described here is not limited to the specific 
surface discretization used here; this discretization was chosen because it has been used in past works on 
the dynamics of fiuid-fiUcd clastic capsules and drops in flow 
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In the present work, the surface of a capsule is discretized into triangular elements. Triangulation of a 
sphere is achieved by mapping the vertices of an icosahedron, which has 12 vertices and 20 triangular faces, 
to the surface of the inscribed sphere [5^ . This procedure will, therefore, give 20 elements on the surface of 
the sphere. Further refinement is obtained by subdividing each triangular face of the icosahedron recursively 
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into 4 equal triangular elements, with all the vertices (and consequently the elements) again being mapped 
to the surface of the inscribed sphere as described above. The number of elements (-/Va) and the number of 
vertices (Nf,) obtained by this procedure can be expressed as Na = 20 • 4*^ and Ni, = Na/2 + 2, where k is 
the level of refinement (fc = corresponds to the original icosahcdron) . Note that the 12 original vertices 
of the icosahcdron have a coordination number 5, while the remaining vertices have a coordination number 
of 6. As an example, a sphere subdivided into A^a = 1280 elements with Ni, = 642 vertices is shown in Fig. 



(4(a) I 



3.3.1. Basis functions over elements 

Linear basis functions are used over each element. All computations over a triangular element is per- 
formed by mapping to or from the parent triangle [26 1. The parent triangle employed in this work is shown 
in Fig. ( |4(b)[ ), where ^ and rj denote the natural coordinates. The basis functions associated with the nodes 
1,2, and 3 are respectively given in natural coordinates by 

0i(e,'7) = l-e-'7, (16a) 

02(e,?7)-e, (16b) 

Mtv)=V- (16c) 
As an example, the position vector x as a function of natural coordinates x(^,?7) is obtained as 

X(e, V) = 4^1 (e, V) Xl + 02 (? , 7?) X2 + 03 (C, ??) X3 , (17) 

where Xi, X2, and X3 are the real space positions of vertices 1,2, and 3 respectively. The same procedure is 
employed to obtain the value of any physical variable (e.g. velocity) at coordinates (^, rf) over the domain 
of the parent triangle. 

3.4. Single layer integral 

Let the single layer integral over the surface S be denoted by 

u.j(z)= / g,(x)Gj,(z,x)d5(x), (18) 
Js 

where 9i(x) is the single layer density, while w(z) is assumed to represent the velocity at point z. In order 
to employ GGEM as discussed in Sec. (j3.2[) to compute the above integral, we write this equation in the 
form 

u;,(z)= / / g,(x)5(y-x)G,,(z,y)d5(x)dy(y), (19) 
J s Jv 

where V represents the volume of the domain and 5 is the three dimensional Dirac delta function. It is easy 
to see that both the expressions for vir(z) in Eqs. ([TS]) and (fT9|) are identical. Next, we write the Dirac-delta 
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function as a sum of the local and the global density introduced in Sec. (|3.2p : see Eq. ([5]). Consequently, 
we have 

w,{z)^ f / g,(x)(p;(y-x)+p<,(y-x)) G,,(z,y)d5(x)dF(y), (20) 
Js Jv 

Next, we separate the integrals associated with the local and global densities, and write the contribution 
due to the local density as 



^.(z)- y^g,(x) (^y^p,(y-x)G,,(z,y)dF(y)j d5(x). (21) 

It is easy to see that the above integral can be written as 

^(z) = / q,(x)G^,(z,x)d5(x), (22) 
Js 

where G' has been defined in Eq. ([T^ . This follows from the fact that the local Green's function G'(z,x) 
can also be constructed by the superposition of Green's function G(z,y) weighted by the density p'(y — x), 
i.e. 

G^,(z,x)=/ p,(y-x)G,,(z,y)dF(y). (23) 
Jv 

It is important to emphasize that the domain was assumed to be unbounded in arriving at Eq. p3[). This 
is always the case for the local problem as discussed in Sec. p.2p : any error in the boundary condition 
introduced due to this assumption will be accounted for in the global calculation. Next, consider the 
contribution from the global density in Eq. ()20[). which we write as 

w'^{z)= JJ^q,{^)p,{y-^)G,,{z,y)dV{y)dS{^). (24) 

It can shown that w^(z) satisfies 

- Vp^"'(z)+/^VV3(z)+/^n^'(z) =0, (25a) 
V-wf(z)=0, (25b) 

where the density 11^ (z) is given by 

nf(z) =87r / q(x)pg(z-x)dS'(x). (26) 
Js 

The boundary condition for the global solution comes from the known local solution (w' (z) for z at the 
domain boundary) and the given overall boundary conditions, such that the sum of the local and the 
global solution satisfies the overall boundary condition; see Sec. p.2p . Having expressed the single layer 
integral in a form suitable for its computation with the GGEM technique, we next describe its numerical 
implementation. This includes the numerical solution of the local problem represented by Eq. (j22p and the 
global problem represented by Eq. 
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3.4-1- Local contribution 

We first consider tlie contribution from ttie local Green's function to the velocity at a given point z, 
which typically is one of the nodes of the elements. We discretize the surface integral in Eq. as 



fc=i •^■^fc 



), (27) 



where the summation is over all the triangular elements N/^Np present in the system, while Sk denotes an 
integral over the element fc. For convenience, all the integrals are performed over the parent triangle. To 
accomplish this, we write the above equation as 

^!(^)= E / / q^i^iLv)) G\,{z,^{£,,i^)) ujd^dr^, (28) 
fe=i -^0 -^0 

where the differential area element dS has been replaced by its equivalent expression 

dS = ujd^dr] = |xj x x^ | c?^ dr]- (29) 

As noted earlier, the value of any quantity can be obtained at coordinates (^, 77) by the usual interpolation 
from the corresponding values at the nodes of the triangle, e.g. see Eq. p7|) . The double integral in Eq. 
(|28p is evaluated using the product of two one-dimensional Gauss-Legendre quadrature rule (one for ^ and 
the other for 7y). This proved competitive in terms of computational cost for a given accuracy with Gaussian 



quadrature rules available for a triangular element 



26[, perhaps due to the fact that the integrands are 



not polynomials. In most cases a 4x4 product rule is found to be sufficient for accurate integration over 
a triangular element. In addition, if the vertex at which the velocity is being computed is a member of 
the triangular element over which the integration is being performed, then the integral in Eq. (j28|) over 



the parent triangle is further transformed to polar coordinates (r, 9) [54| . This transformation makes the 
integrand non-singular and hence ensures sufficient accuracy with the same low order product integration 
rule discussed above. Lastly, we note that for computing the contribution to velocity at any given point 
z due to the local Green's function, only triangular elements within a distance of rcut ^ ct^^ from the 
point z need to be considered (typically, rcut = 4/a). This is justified due to the exponentially decaying 
contribution from the integral over an element at separations larger than 0{a^^) from the point of interest. 
As mentioned in Sec. p.2.ip . the near neighbor list required for the local calculation is generated in 0{N) 
time via the cell- linked list method 
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3.4-2. Global Contribution 

Our goal here is to find 11^(2) in Eq. (P5|) . for wliicfi we need to compute tlie integral in Eq. (pS)) . We 
being by discretizing the integral in Eq. (|26[) and write it as 

nf(z)= ^ / q,(x)p3(z-x)d5(x). (30) 

fc=l 

The above integral can be evaluated with any desired quadrature rule, though, due to the smoothly varying 
nature of the integrand a simple trapezoidal rule proves sufficient. Note that the trapezoidal integration 
rule essentially reassigns the contribution from the surface of the triangular clement to its three vertices in 
equal proportions. Consider first the integral over an element Sk in Eq. pop . which as per the trapezoidal 
rule is expressed as a sum of contributions from its three vertices as 

[nf(-)]..-E(^^^^) P.(z~x^>), (31) 

where Ak is the area of the triangular element Sk, p denotes the vertex number of the given element fc, x'^p 
denotes the coordinate of the p*'' vertex of the triangular clement fc, and [11^ (z)]^^ denotes the density at z 
due to the integration over the element Sk only. The term in the parenthesis in Eq. pi[) can be considered as 
the strength of the global density at the node x'"'p due to the clement fc; by summing it over all the elements 
k to which a given node belongs (let us say that this node is globally represented by x^), one obtains the 
total strength of the global density at this node, say q'' . The overall density at a point z is then obtained 
by adding contributions from all the nodes present in the system, 

nf(z)=E<zf'p,(z-x^), (32) 

where iVj, is the total number of nodes in the system. We also note that due to the exponentially decaying 
nature of Pg{z — x'') as a function of distance from x**, we consider only those nodes for computing the 
density at a point z which arc within a distance rcut ^ from it. Once n^(z) is evaluated, we solve 
the set of equations in ([^5)) using the procedure described in detail in Sec. p.2[) . This gives us w3(z) at 
the mesh points. The velocity w^(z) at any point not on the mesh is obtained using 4*'' order Lagrange 
interpolation. Once w^{z) is known, the overall single layer integral w(z) is obtained as: 

w(z) = w'(z) +w»(z). (33) 
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3.5. Double layer integral 

We now describe the evaluation of the double layer integral. We denote the double layer integral over a 
surface S (similar to that in Eq. ^ by 

z;,(z)=nfc(z) / g,(x)T,,fc(z,x)d5(x). (34) 
Js 

If we define a stress tensor cr(z) as (note that one needs to multiply it by /i to get units of stress) 

crjfc(z) = / (7j(x)Tjifc(z,x)dS'(x), (35) 
Js 

then we have the following relationship between v(z) and o'(z): 

?;j(z) = nfc(z)a-j7,(z). (36) 

The motivation for introducing the stress field cr(z) in Eq. ([55]) should be clear now, as it is the stress field 
associated with the following velocity field 

Wj{2,) = / g,(x) G,,(z,x) dS{-x.). (37) 
Js 

As described in the previous section p.4p . we write the above velocity field as the sum of a local w'(z) 
and global w3(z) velocity fields, and denote the corresponding stress fields by cr'(z) and erf (z) respectively. 
Using (j36|) . we obtain the corresponding local and global contributions to the double layer integral as 

i/.(z) ==nfc(z)<7j,(z), (38a) 

^;|(z)=n,(z)af,(z). (38b) 

We describe the computation of the global contribution to the double layer integral first, as it is a straight- 
forward extension of the procedure presented in Sec. p.4p . This will be followed by a discussion of the 
procedure for computing the local contribution to the double layer integral. 

3.5.1. Global contribution 

Here we describe the procedure to compute the global contribution to the double layer integral. Consider 
the global velocity field w^(z) and the pressure field ^i^" {z.) associated with the global force density 11^ (z); 
see Eq. (j25|) . The procedure to compute the velocity and pressure field for this global distribution of density 
has been discussed in detail in Sec. p.4[) . Once these arc known, one can obtain the stress field cr^{z) from 
the usual Newtonian constitutive equation as 
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As was mentioned in Sec. p.2p . the differentiations required in the above expression are performed in the 
transform space 0,14^; consequently, the stress field is known with spectral accuracy at the mesh points. 
Once the stress tensor is obtained at the mesh points, the stress at the nodes of the elements are obtained 
using 4*'' order Lagrange interpolation. Lastly, (z) is obtained from Eq. (j38bp . 

3.5.2. Local contribution 

Consider the velocity field due to the local density, which is written as 



w\^)^ / g,(x)G^,(z,x)d5(x). (40) 
Js 

It is trivial to show that the stress field associated with the above velocity field is given by the following 

ajfe(z)= / g,(x)rj,fe(z,x)d^(x), (41) 

and consequently the local contribution to the double layer integral is given by 

v\{2) = n,(z) / g,(x)Tj,,(z,x)d5(x). (42) 
Js 

The above integral is discretized in the same fashion as in Sec. p.4p . and is written as a sum of integrals over 
the triangular elements. Again, as in Sec. p.4.ip . all the integrals are performed over the parent triangle 
as follows 



v]iz) = nk{z) / q^ i^iC V)) Tj,,iz,^{tv)) ^d^drj. (43) 

Jo Jo 



3 

k=l 



The integral in the above equation is performed with the product of two one dimensional Gauss-Legcndre 
quadrature rule; see Sec (|3.4.ip for details. In addition, when the vertex at z is a member of the element 
over which integration is performed, the integrand is singular and requires special treatment. At this point 
it is worth noting that the double layer integrand has the same 1/f singularity as the single layer integral 



47 1, as X • n(xo) ~ for small f. So in principle one may employ the polar coordinate transformation to 



regularize the singular double layer integral. In practice, because we employ flat elements, the condition 
X • n(xo) ~ is not valid. Note that the normal vector at any given vertex is taken as the area averaged 
normal vector of triangles to which it belongs, from which it follows that the numerical double layer has 
a singularity. This stronger singularity can be avoided by employing surface elements that yields a 
continuously varying normal vector such as splines {sgI ]. 

In the literature, desingularization of the double layer integral is usually achieved by singularity sub- 
traction [sot . This desingularization scheme is also applicable for the current formulation p4p . though this 
procedure will require 9 times the computational effort of performing the double layer integral itself; we omit 
the details here. We therefore look elsewhere for more computationally efficient schemes for performing the 
singular part of the double layer integral. The first point to note is that the singularity in the stress tensor 
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T(z,x) is contained entirely in its local part T'(z,x). Moreover, desingularization is necessary only when 
the target point z is one of the vertices of the element over which integration is being performed. If this is 
the case, we proceed by replacing the normal vector n(z) outside the integral in Eq. (j43|) by a vector n(x) 
inside the integral, which leads to the following expression for the integral over the current element k: 

[v'^iz)]^^= [ g,(x)Tj,,(z,xK(x)rf5(x). (44) 

The vector n(x) at any point x(^, 77) on the element is defined as per the following equation: 

n(x) =nA(/)i(e,??) + n(z),^2(e,^)+n(z),^3(e,?7), (45) 

where ha refers to the normal vector of the current triangular element. In writing the above equation, we 
have assumed that the target point z is mapped to the vertex labeled 1 of the parent triangle; see Fig. |4(b)[ 
It is immediately clear from Eq. (|45)) that n(x) will take the value ha when the source point coincides 
with the target point (i.e., when x = z), while it will tend to n(z) as the source point x moves away from 
the target point z over the current element. Now, if we substitute the expression of T' from Eq. (|12l) into 
Eq. ([33]), one immediately sees that the singular parts of the integrand tends to zero as the source point 
approaches the target point since x • ha = 0. This is due to the fact that elements are flat and x lies in the 
plane of the element, while ha is normal to the element. In order to estimate the error introduced due to this 
approximation, one can show by simple Taylor series expansion that ||n(x) — n(z)|| ~ ||Vn(z)|| hs, where 

— 1/2 

is the characteristic surface mesh spacing scaling as . In addition, since this approximation is applied 

only when the target point is a member of the triangle over which the integration is being performed, the 
previous error gets multiplied by the area of the triangle which is 0{N^^). Therefore, we estimate the error 
introduced in the solution due to this approximation as 0{N^^^^). This completes the evaluation of the 
local contribution to the double layer integral. The total double layer integral is then obtained as the sum 
of the local and the global parts as: 

v,{z)=vUz)+v^^{z). (46) 



4. Membrane Mechanics: Hydro dynamic traction jump 

The solution of the boundary integral equation (|3]) requires the knowledge of the hydrodynamic traction 
jump across the interface Af . This jump in traction is obtained from the membrane equilibrium condition 
as discussed next. Consider first a patch of element on the membrane's surface as shown in Fig. ([S^). 
The forces acting on this patch are the hydrodynamic stresses on the inner and the outer surface and the 
membrane tension at the boundary denoted by contour C. Now, let the membrane tension tensor be given 
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(a) A patch on the membrane (b) A patch on the discretized 

surface 

Figure 5: Defining the contour C enclosing an area (hatched) containing the point of interest P. b is the in-plane normal to 
the contour C . 

by T, then the force balance on the membrane patch is given by 



MdS+ / b • T = 0, (47) 

Sc Jc 

where Sc denotes the area enclosed by the contour C. Using the divergence theorem, one can convert to 
contour integral to a surface integral, which in the limit of infinitesimal area yields 

Af = -V, • T, (48) 



where is the surface divergence operator |3|. The above equation (j48[) has been directly employed in 
several boundary integral implementations to obtain the traction jump Af ; see e.g. Lac et al. [sSj. For flat 
elements as in this work, we note that t is constant over each clement such that its surface divergence is 
identically zero, though there is a jump in its value across elements and consequently the contour integral 
in the Eq. (|T7)) is generally expected to be non-zero. The Contour integral can therefore be used to obtain 
the traction jump as js^ 

Af / h-Tdl, (49) 

where the Ac is the area enclosed by the contour C; Fig. ^p) shows an example of the contour for the 
discretized surface. To proceed with our implementation, we first interpret the contour integral on the right 
of (j49p as the reaction force on a given node obtained under the condition that the entire elastic energy 
stored in the particle membrane has been reassigned to the vertices of the discretized triangular elements. 



y 



to compute the reaction force 



We then use the principal of virtual work as presented by Charrier ct al. 
at the vertices. Once the total reaction force Fp at a given vertex P is known, we obtain the traction 
discontinuity at that vertex as 

Afp = -^ (50) 
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where Ap is area assigned to the vertex, which is essentially the area enclosed by a hypothetical contour 
around the vertex P; see Fig. ([SJd). We call the contour hypothetical as we never explicitly define it here. 
For the area assignment to the vertex, we use a very simple rule where each vertex of the triangular element 
is assigned a third of the triangular element's area. Therefore the total area Ap is obtained as (1/3) of 
the total area of the triangular elements of which the given vertex P is a member. We believe the method 
outlined here is substantially simpler to implement than employing the contour integral explicitly. In the 
remainder of this section, we outline the procedure employed for computing the reaction force at the vertices 
of the triangular elements. 

We begin by introducing the formalism for describing the kinematics of the membrane deformation. This 
formalism is mostly clearly presented for deformations in a plane, which for the moment is taken to be the 
xy plane. Let (x, y) and {X, Y) denote respectively the undcformed and deformed coordinates of a material 
point, with respect to a fixed set of Cartesian axes. If u and v denote the displacements of the material 
point in x and y directions respectively, then 



X = X + u, 
Y = y + V. 



(51) 



The relationship between an infinitesimal line segment before and after the deformation can be expressed as 

(52) 



dX 




1 + du/dx 


du/dy 




dx 


dY 




dv/dx 


1 + dv jdy 




dy 



or compactly as 



dX = F • dx. 



(53) 



(54) 



where F is the deformation gradient tensor. The square of the distance between the two neighboring points 
after deformation is given by 

dS^ = dX • dX = dx ■ G ■ dx, 
G = F^ F, 

where G is a symmetric positive definite matrix. We denote the eigenvalues of the G by and A^, such 
that Ai and A2 are the principal stretch ratios. For a thin membrane that displays no resistance to bending, 
the strain energy density W of the membrane is a function of Ai and A2. Here we consider the capsule to 
be an infinitely thin neo-Hookean membrane, for which the strain energy density is defined as P| 

1 



G 
2" 



A^ + A^ 



(55) 



Here G is the two-dimensional shear modulus for the membrane, having units of force per unit length. To 
compute the reaction force at the nodes, we adopt the finite element approach of Charrier et al. [8|. Only 
the briefest account will be given here; for details the reader is referred to the original reference. In the 
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approach of Charricr et al. [8[ , the membrane is discrctized into flat triangular elements such that the strain 
is uniform over an element. Moreover, it is assumed that an element remains flat even after deformation. 
The forces at the nodes are then determined from the knowledge of the displacement of the vertices of the 
element with respect to the undcformed element followed by the application of the principal of virtual work, 
such that the computed forces and the known displacements are consistent with the strain energy stored 
in the element. For an arbitrarily oriented element, rigid body rotations and translations can be defined 
to make the deformed and undeformed state in the same plane. Note that the rigid body rotations and 
translations have no effect on the strain energy and consequently the forces. The forces are then computed 
using the coplanar formalism discussed above. Finally, these forces are transformed back to the frame of 
reference in the deformed state by applying the inverse transformation. The total reaction force at a node 
is obtained as a sum of reaction forces at that node due to contributions from all the triangular elements 
of which it is a member. Once the reaction force at any given node is known, the hydrodynamic traction 
discontinuity at that node is obtained from Eq. ([50)) as detailed earlier. 

5. Solution procedure and parameters 

In this section we describe solution methods and parameters for the equations resulting from the formu- 
lation just presented. 

The heart of the computation is the determination of the fluid velocity at the element nodes using 
equations ^ and ([3]). To compute the velocity from Eq. ([2|), we first need to compute the single layer 
density q(x). This is obtained from the solution of the integral equation ([3]). Upon discretization of the 
double layer integral in Eq. ([3]), which was discussed in Sec. p.Sp . we obtain a linear coupled system of 
equations for q^, where is a vector of length 3Ni, denoting the value of q(x) at the element nodes; Nf, is 
the number of clement nodes in the system. We express this linear system of equations as 

(I+£.D^).q'' = b, (56) 

where k is a diagonal matrix of size 3Nb x 3Nb denoting the value of Km in Eq. ([3]) at each element node, 
denotes the discrctized double layer operator of size 3 A'';, x 3Ni,, while b is a 3Ni, vector denoting the 
known right hand side of Eq. at the element nodes. The above system of equations is solved itcratively 
using the GMRES algorithm [ssj. An important benefit of this iterative procedure is that the matrix in 
the parenthesis above is never explicitly computed; at each iteration step only the product of the above 
matrix with a known vector generated by the algorithm is to be computed. The procedure to compute this 
matrix vector product is similar to computing v(z) in Sec. p.Sp (see Eq. I34p at the element nodes for a 
known q*". We take the initial guess for q'' cither from the previous time step, or from the previous stage 
if a multistage method is employed as is the case here. This leads to a substantial savings in the number 
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of iterations required for convergence. Iterations were terminated when the L2 norm of the current residual 
vector S relative to the norm of the right hand side b was less than 10~^; the residual vector S is defined as 

S = (I+£.D^).q^-b. (57) 

In addition to the convergence in the residual vector S, which is obtained naturally as part of the iterative 
procedure, it is also important to investigate the convergence in the solution; the error in solution denoted 
by vector Sg is defined as 

S, - q" - qL (58) 

where we have defined the "exact" solution as q^^, which can be obtained by solving the above system of 
equations to a very small tolerance, e.g. 10"^". The above tolerance of 10"'^ for the residual vector relative 
to the norm of the right hand side leads to an error of the same order for the error vector relative to the 
exact solution for well conditioned systems. This implies an accuracy of 0.01% for most cases. Even for 
the worst cases in the present work, the relative error in the solution was always less than 0.1%. In most 
cases convergence was achieved in less than 5 iterations; see Sec. ([7]) for some examples. For systems with 
high viscosity contrast and/or with large number of particles, the matrix may become ill-conditioned and 
a preconditioner may become necessary. In the present study, no prcconditioncr was employed, though 
multigrid preconditioners for Stokes fiow may be useful [l^ IsGj . 

The iterative procedure described above gives the q'' vector at the element nodes. These are subsequently 
substituted in Eq. ^ to compute the corresponding velocity at the element nodes; the numerical procedure 
described in Sec. p.4p is employed to compute this single layer integral. We denote the velocity thus 
computed at the element nodes by a iN^ vector u**, which is used to evolve the position of element nodes 
x'' as per the equation 

u''. (59) 



dt 

The time integration in the above equation is performed via the second order midpoint method, which 
belongs to the family of explicit Runge-Kutta integrators [s^l ■ The time step At employed in this work was 



set adaptively using the rule 
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531: 

7At = 0.5 Ca/i^/a, (60) 



where /i™ is the minimum node-to-node separation in the system (which does not have to be for two points on 
the same particle), a quantifies the length scale of the particle such as the radius for a spherical particle, 7 is 
the shear rate, while Ca is the capillary number, expressing the ratio of viscous and interfacial stresses (for a 
spherical capsule with radius a and shear modulus G, the capillary number is defined as Ca = fi-ya/G). This 
rule gave a stable evolution with time; a time step twice this value led to an instability in some simulations 
presented later, while a time step half this size gave nearly indistinguishable result. The volume of the 
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Figure 6: Schematic of the slit geometry for the test problems in Sees. 116. Il l and 1 16. 2| I. A single sphere, either a rigid particle 
or a drop, is placed in a slit geometry with the channel height being H. The radius of the sphere is a, while its center is at j/c- 
A pressure driven flow is considered with Uq being the center line velocity. The confinement ratio is 2a/ H. Periodic boundary 
conditions are employed in x and z directions, with periodicity being and Lz respectively. 



particle was found to be well conserved with time. For example, for a capsule with A = 5 in shear flow at 
Ca = 0.6, the volume changed by an average of 10^^ of its original value over a unit strain. In the present 
work, a volume correction was performed only when the volume of the particle deviated by more than 10"'^ 
of its original value; the procedure employed for this correction is described in 

Finally, we consider parameters related to the numerical solution procedure. Many of these have already 
been introduced earlier, and are repeated here for completeness. The first parameter is TVa, which refers 
to the number of triangular elements employed to discretize the surface of each of the Np particles in the 
system; the number of element nodes per particle is denoted by Nf,. Next, there are several parameters 
associated with the GGEM methodology described in Sec. (|3.2[) . First is the length scale associated 
with the quasi-Gaussian global density. The local solution as well as the assignment of the global force 
density to mesh points is truncated beyond a distance of Tcut = 4/a from the origin of the singularity. The 
solution of the global problem (Sec. 13. 2p requires one to define a three dimensional mesh with N^, Ny, and 
N2 mesh points in x, y and z directions respectively. This gives a mean mesh spacing in the three directions 
as — Lx/Nx, A.ym = H/{Ny — 1) and Az„i = Lz/N^] all three mean mesh spacings are kept equal 

unless otherwise mentioned. An important parameter denoting the resolution of the GGEM methodology 
is aAy,„, with the resolution and hence the accuracy of the method increasing with decreasing aA?/„j. As a 
rule of thumb, we require aAy„i < 1; see Sec. p.2.3p for details. Choices of above parameters are specified 
below in the descriptions of the test problems. 
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(c) Effect of the spatial period L 

Figure 7: Rigid particle in a slit: (a) Comparison of the translational velocity of a rigid sphere in a slit geometry with results of 
Staben et al. [5^ . The confinement ratio of the particle was 2a/ H = 0.6. The horizontal axis gives the velocity of the particle 
non-dimensionalized by the centerline velocity of the undisturbed fluid, while the vertical axis gives the location of the center 
of the sphere along the gradient direction non-dimensionalized by the height of the channel. Walls are present at y = and 



y = H. Simulation parameters were 
for details. 



: 61, aAy„ 



0.44, A^A = 5120, Lx = Lz = 5H, and Axm = Aym = Azml see text 
(b) Convergence of the velocity for yc/H = 0.5 with oAj/m- For this calculation Vcut = 0.5a is kept fixed, while 
(c) Effect of the spatial period L = Lx = Lz on the translational velocity of a particle with 
yc/H = 0.5. Also shown is a linear fit to the data in the plot. All simulation parameters in (c) were the same as in (a) above 
except for Lx and Lz which were varied. 



6. Numerical Results: Validation 

6.1. Single layer validation: Rigid particle in a slit 

As a validation of the single layer integral implementation, we consider a rigid sphere between two parallel 
walls and subject it to a pressure driven flow with a centerline velocity C/q as shown in Fi g. (l6l) . For a rigid 
particle, the velocity at a point xq on the surface satisfies the following integral equation |3C| 



wj(xo) = u°°(xo) - / /i(x)G-,i(xo,x)(i5'(x), 



(61) 
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where f is the traction on the external surface of the sphere due to the stresses in the fluid. The velocity at 
the surface of the particle u(xo) in the above equation can be written as 



u(xo) = U + n X (xo -Xc), (62) 

where U and $1 represent the translational and rotational velocities of the particle, and Xc denotes the center 
of the sphere. Here we take the particle to be force and torque free and our goal is to compute the velocity 
and angular velocity of the particle. For this particular problem, the surface of the sphere was discretized 
into A^A = 5120 triangular elements with Nf, = 2562 vertices. The unknowns in the discretized system are 
SNh tractions at element vertices, along with U and H. The force and torque free condition along with the 
discretization of Eq. (|5T|) gives 3Nh + 6 equations, which were solved iteratively using the GMRES algorithm 

Translational and rotational velocities of a rigid sphere between two infinite parallel walls have been 
reported previously by Staben et al. [58[. Here, we compare the translational velocity obtained in the 
present work with their results for a fixed confinement ratio of 2a/ H = 0.6 and for various positions of the 
sphere's center along the channel height j/c- This comparison is shown in Fig. ([7^), where the velocity of the 
particle has been non-dimensionalized by the velocity of the undisturbed flow at the centerline Uq, while the 
height of the sphere's center has been non-dimensionalized by the channel height H. Very good agreement 
between the two results was observed, with the discrepancy typically being less than 0.8%; the source of 
the slight discrepancy is discussed below. The GGEM parameters employed in the above calculation were: 
Ny = 61 and Vcut = 0.5a, which gives aAy™ = 0.44; note that — A/a. The mean mesh spacing was 
equal in all three directions (Ax™ = ^y-m = Az^), and the spatial period in both x and z directions were set 
to five times the wall spacing: = ~ 5H. Convergence of the particle velocity with respect to aAj/m is 
demonstrated next in Fig. ([TId) for a particle placed at the centerline, i.e. yc/H — 0.5. For this calculation 
fcut = 0.5a was held constant, while aAy™ was varied by varying Ny between 17 and 81 {N^ and varied 
between 80 and 340). As could be seen in the figure, the velocity of particle reaches its converged value for 
a Ay™ < 0.5 and shows very little variation with any further increase in mesh resolution. Also shown in 
this plot is the result of Staben et al. 58| which reveals that the velocity obtained in this work converges 
to a slightly lower value than the previous reference. The source of this discrepancy can be traced to the 
periodic boundary conditions employed in x and z directions in the present work; Staben et al. 



used an 



unbounded domain in these directions. In the present case, we can easily estimate the result for an infinite 
box by observing its trend in a series of simulations with varying spatial period L (L = ~ L^). This 
procedure is commonly used in triply periodic simulations to remove the effects of the periodic boundary 
conditions; see, e.g. [3j|. In this particular example, we numerically find that the periodic image effects 
decay as as shown in Fig. where we have plotted the translational velocity of the particle against 
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Figure 8: Volume averaged translational velocity of a spherical drop in a slit. The confinement ratio of the drop is 2a/ H = 0.6, 
while its center is at yc/H = 0.4. The simulation parameters are: Ny = 61, aAym = 0.44, TVa = 5120, and Axm = Aj/m = 
A^m. (a) Effect of the spatial period L = Lx = on the drop's velocity for different viscosity ratios A. Note the velocity has 
been non-dimensionalizcd by Uac, which is the estimated velocity of the drop as L — > oo. This is obtained by fitting a straight 
line to the data which have not been non-dimensionalized by Uao', the y-intercept of this fit gives Uoo- (b) Comparison of the 
volume averaged translational velocity of a drop as a function of viscosity ratio with the results of Janssen and Anderson [28| 
(JA). The velocity in the current work corresponds to Uaa in (a) above. 



L^^. The y- intercept of the hnear fit through the data points in the previous plot then gives an estimate 
of the particle velocity in an infinite slit. This value comes out to be 0.871C/o (rounded to three significant 
digits), which is exactly equal to the value reported by Staben et al. [58 1. Recall that their boundary integral 



formulation is based on the slit Green's function of Liron and Mochon 
with those of Staben et al. 



42|; the agreement of our results 



thus implicitly validates our Green's function implementation for a slit with 



421 . In addition to the periodic image effects, 



respect to the Green's function provided by Liron and Mochon 
a slight discrepancy between our results and those of Staben et al. [ssj can also be expected in cases where 
the particle-wall separation is very small. In such problems, a large lubrication pressure can develop in the 



region around the small gap 



3C 



33| . To obtain accurate solutions in this case, the surface discretization of 

. No fundamental 



the particle near the small gap must be adaptively refined as was done by Staben et al 
changes to the present formalism would be required to implement adaptive refinement. 

6.2. Single and double layer validation: A Drop in a slit 

Having validated the single layer integral, we next move on to the validation of the double layer integral. 
For this, we consider the same geometry and bulk fiow as in the above test case (Fig. [B]), but now consider 
a spherical drop instead of a rigid sphere. The motion of the drop can be obtained by first solving Eq. ([3]) 
for q(x), which upon substitution in Eq. ([2]) gives the velocity on the surface of the drop. We first point out 
that for a spherical drop, the interfacial traction jump Af is inconsequential. This is due to the fact that 
Af is uniform in strength and acts radially everywhere, which when combined with the incompressibility of 
the fiuid implies zero velocity contribution from this term. Once the velocity at the surface of the drop is 
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known, we compute the volume averaged velocity of the drop as 

Ui^ jz ( u^dV = I {ujnj)xidS, (63) 
V Jv ^ J s 

where V represents the volume of the drop, while n is the unit normal vector at the surface. We computed the 
instantaneous volume averaged velocity of a spherical drop placed at yd H ~ 0.4 and with a confinement ratio 
of 2a/ H = 0.6 for different drop viscosity ratios A. The simulation parameters for this calculation were kept 
the same as in the previous section, namely Ny = 61, TVa = 5120, aA?/,„ = 0.44, and Aa;,„ = A?/,„ = Az^. 
Just like in the case of rigid particle above, we again find the same L^^ scaling of the periodic image effects 
on the drop velocity. This is shown in Fig. ([5^) for drops of different viscosity ratios A, where the velocity 
of the drop has been non-dimensionalized by the corresponding velocity estimated for an infinite slit using 
the procedure outlined above in Sec. (j6.ip . Interestingly, for a drop with A = 1, there is no observable 
periodicity effect, while the drop velocity increases with increasing L for A > 1 and decreases with increasing 
i for A < 1. In Fig. (|5Jd), we compare the results in the present work corrected for periodicity effects 
with those of Janssen and Anderson [28[ for several different viscosity ratios (A S [0.5, 1, 2, 5]). A very good 
agreement between our results and those reported by Janssen and Anderson 



28l | is evident at all viscosity 



ratios (error < 0.15%), thereby validating our implementation of the single and the double layer integral. 

6.3. Validation of the capsule membrane mechanics and the overall implementation 

We next consider a capsule in a simple shear flow. To enable comparison with literature results in 
an unbounded domain, we consider a large simulation box with ^ H = = 15a, where a is the 
radius of the initially spherical capsule placed at the center of the box. Other simulation parameters were: 
Ny — 97, rcut — 1-0, aAym = 0.625, A^a = 5120, and Aa;j„ = A?/™ = Az,„; a convergence study with 
these parameters will be presented later in this section. The system is subjected to simple shear flow and 
we follow the evolution of the shape of capsule at various capillary numbers (Ca = pL'^a/G) and viscosity 
ratios (A). Lac et al. [35| showed that a membrane lacking bending resistance buckles at high or low Ca; 
the origin of this buckling has been shown to be numerical |39j . We saw a similar behavior and therefore 



restrict the results reported here to 0.3 < Ca < 0.6. In this regime, the capsule shape evolution appeared to 
be stable with no apparent buckling. To characterize the shape of the deformed capsule, we introduce the 
commonly employed Taylor deformation parameter D defined as 

where L and B are the maximum and the minimum distance in the shear plane of a point on the surface of 
the capsule from its center. We use this as the definition of D, though some authors, e.g. Ramanujan and 
Pozrikidis 5J], instead find a triaxial ellipsoid with the same inertia tensor as the given capsule, and then 



take L and B as the major and minor axis of that ellipsoid. 
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Figure 9: A = 1: Time evolution of the Taylor deformation parameter D at (a) Ca = 0.3 and (b) Ca = 0.6. Lac et al in the plot 
refers to the results of Lac et al. [35l . while RP refers to the results of Ramanujan and Pozrikidis [541 . Simulation parameters 
were: Lx = Ly = H = 15a, Ny = 97, rcut = a, aAym = 0.625, A'^a = 5120, and Axm = Aym = Azm- (c) Convergence of 
D at Ca = 0.6 with respect to aAym- In this study rcut = a was held fixed, while Ny was varied, (d) Convergence of D at 
Ca = 0.6 with respect to A^a- Rest of the parameters are the same as in (a) and (b). 



The evolution of the deformation parameter D for a capsule with unit viscosity ratio (A = 1) is shown 
in Figs. and (j^b) at two different capillary numbers Ca. The time in the figures has been non- 



Mi 



dimensionalize d by the shear rate. For comparison, D values reported by Ramanujan and Pozrikidis 
and Lac et al. [ssj are also plotted. Note that Ramanujan and Pozrikidis [s^ have used a zero-thickness 
shell model for their capsules, though that gives only marginally lower deformation than a neo-Hookean 
capsule at the same Ca [sJl- Moreover, they used the Young's modulus for computing their Ca, such that 
our results should be compared with their results at a Ca which is (1/3) of the Ca in this work. The data 
in Figs. ([9^) and both show that the evolution of D in this work is in very good agreement with the 
correspon ding results of Lac et al. [35| . Our results are also close to the values reported by Ramanujan and 



Pozrikidis 



54| . though the latter consistently display a slightly lower D. The broad agreement of D between 



our values and the literature values validates our implementation of the membrane mechanics along with 
other aspects of our method such as time-stepping and the already validated single layer integral. 
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Figure 10: A = 5: Time evolution of the Taylor deformation parameter D at (a) Ca = 0.3 and (b) Ca = 0.6. Results from 
the present study are compared with the work of Ramanujan and Pozrikidis 1541 (RP) and Le and Tan [ssll (LT). Simulation 
parameters were the same as in Fig. (|9]l. (c) Convergence of D with A^a at Ca = 0.6. 



Next, we demonstrate the convergence of the steady state D at Ca = 0.6 with respect to the GGEM 
parameter aA?/,„ in Fig. (jHt)- For this calculation, rcut = o, was held fixed, while Ny was varied between 65 
and f f3. As could be seen, a convergence in D is observed at Ny ~ 97 corresponding to aAy™ = 0.625 - 
the parameter set for which all results are presented in this section. The convergence of the steady state D 
with respect to the number of triangular elements A^a is demonstrated in Fig. (jOji) for a Ca ~ 0.6 capsule. 
The three data points in this plot correspond to simulations with Aa = 320, 1280 and 5120 elements. It is 
clear from this plot that the solution converges linearly with A^a^' which is expected for linear elements. 

We now briefly discuss the deformation parameter results for capsules with non-unit viscosity ratios 
(i.e. A 7^ 1). These are reported in Fig. (fTOj) for A = 5 and in Fig^^ (fTTj) for A = 0.2, and are compared 



with the results of Ramanujan and Pozrikidis j54l | and Le and Tan [38j, with the latter reference employing 
the immersed boundary technique for their simulations. At A = 5, the results for D in the present work 
are slightly smaller than those reported by the previous authors at Ca = 0.3. At Ca = 0.6, an excellent 
agreement with the results of Ramanujan and Pozrikidis 
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54| is observed, though our results for D are slightly 
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Figure 11: A = 0.2: Time evolution of the Taylor deformation parameter at Ca = 0.3 and (b) Ca = 0.6. Results from the 
current study arc compared with the work of Ramanujan and Pozrikidis [sSl (RP) and Lc and Tan |38|| (LT). Simulation 
parameters were the same as in Fig. (|9]l. (c) Convergence of D at Ca = 0.6 with A^a- 



lower than those of Le and Tan 



38|. We next show the convergence of the steady state D at Ca = 0.6 with 



respect to in Fig. (fTUb ). These resuhs indicate that the error decays as N^^^'^, thereby implying that 
the error incurred in the calculation of the singular double layer integral dominates the overall error; see 
Sec. p.5.2p for details. 

We next discuss the results for the deformation parameter for capsules with A = 0.2 (Figs. [TTk andfTTb). 



381 at both Ca ~ 0.3 and 



In this case, wc observe a very good agreement with the results of Le and Tan 

Ca = 0.6. The values reported by Ramanujan and Pozrikidis [3] are also close, though they are marginally 
lower than the results in the present study. Lastly, we show the convergence of the steady state D at 
Ca = 0.6 with respect to Na in Fig. (fTTb '). In this case the error is observed to decay at the expected rate 
of - this probably implies that the error incurred in the calculation of the singular double layer integral 
is dominant only at high values of A. 
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Figure 12: Pair collision: A = 1. (a) Shows the separation between the centers of mass of capsules in the gradient direction as 
a function of their separation in the flow direction. Also plotted is the corresponding result from Lac et al. [3^ . Simulation 
parameters were: Lx = Ly = H = 30a, Ny = 129, rcut = 2a, aApm = 0.469, A'^a = 1280, and AcCm = ^y-m = Azm. 
(b) Convergence of the maximum displacement of either particle (absolute value) from its initial position along the gradient 
direction {y) as a function of aAym- Data points in this curve were obtained by holding rcut = 2a fixed and varying Ny. (c) 
Convergence of the maximum absolute displacement of a particle in the gradient direction with N/\ . 



6.4- Pair collision 

As a final test problem, we consider the collision between a pair of capsules with A = 1 in a simple shear 
flow and compare the results with the work of Lac et al. [36| in Fig. p^ ). We first describe the problem 



setup. The size of the cubic box (slit) for this problem was set to 30a to approximate an unbounded domain; 
simulations were performed with Ny ~ 129, rcut ~ 2a, aAy^ = 0.469, TVa = 1280, and Ax,„ = A?/,„ = Azm- 
The two capsules were initially kept in the same flow-gradient plane [x ~ y), such that the initial separation 
in the flow direction was X2 — xi — —8a, while the initial offset in the gradient direction was y2~ Hi — 0.5a. 
As in Lac et al. [s^, we preinflate the capsule by 5%, i.e. the radius of the spherical capsule was increased 
by 5% over its rest value, and this new increased radius is denoted by a. For a spherical shape, this inflation 
does not lead to any flow due to the incomprcssibility condition as discussed above in the case of a spherical 
drop. Nonetheless, this inflation keeps the membrane in a state of tension at rest and, if sufficient, will 
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(a) Snapshot for A = 1 (b) Snapshot for A = 5 




(c) Convergence (d) A dependence 



Figure 13: Suspensions of capsules at Ca = 0.5 and volume fraction ip = 0.15: (a) Snapshot of A = 1 capsule suspension with 
Np = 120, (b) Snapshot of A = 5 capsule suspension with Np = 120, (c) Convergence of the apparent intrinsic viscosity for 
A = 1 capsules in Np = 120 particle system. The simulation parameter aAym was 0.5 for the run labeled 120, while it was 
0.417 for the run labeled 120* , (d) Effect of A on the apparent intrinsic viscosity for A^p = 120 particle system. The simulation 
parameters in all cases above were: Vcut = o,, aAym = 0.5, A^a = 320, and Axm = Aym = Azm unless otherwise mentioned. 



prevent buckling during the course of the cohision. With these prehminaries, we return to Fig. (jl2b ') where 
we show the relative separation between the center of masses in the gradient direction as a function of the 
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corresponding separation in the flow direction. A very good agreement with the results of Lac ct al 
is evident. Next, we show the convergence of the numerical scheme with aA?/„i in Fig. (jl2b ). where we 
plot the absolute value of the maximum displacement in the gradient direction for cither particle. For this 
calculation Vcut = 2a was held fixed, while Ny was varied between 97 and 145. A convergence is seen at 
aAym = 0.469 corresponding to Ny = 129 ~ the parameters for which results are reported in this section. 
The convergence of the maximum displacement in the gradient direction with respect to is demonstrated 
in Fig. (fT2b '). which confirms the expected error decay rate of N^^. 

7. Multiparticle Simulations: suspension apparent viscosity and computational complexity 



In this section, we report results from several large scale simulations on multiparticle suspensions of 
capsules. In the first part of the section, we present results for the suspension viscosity and discuss its 
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dependence on the viscosity ratio. In the second part, we discuss the expected computational complexity of 
the algorithm and verify it with timing results from the multiparticle simulations. 

7.1. Suspension viscosity 

We consider here a suspension of Neo-Hookean capsules in a cubic slit (Fig. [T3k.b) . These suspensions are 
subjected to a simple shear flow at a capillary number of Ca = 0.5; the volume fraction of the suspension 
is = 0.15, which is typical of the blood flow in the microcirculation [l7|. Capsules with five different 
viscosity ratios A arc considered: A = 1, 2, 3, 4, and 5. The surface of each of the capsules was discrctized 
into A^A = 320 triangular elements, while four different system sizes were considered with the number of 
particles being Np = 15, Np = 30, Np = 60, and Np = 120. In each of the problems, Vcut = a was kept fixed, 
where a is radius of the spherical capsule at rest. Similar to Sec. (|6.4p . the capsules were preinflated by 5% to 
prevent membrane buckling; the radius after preinflation is denoted by a. The number of mesh points Ny for 
the global solution varied between 61 and 121 such that Q!A?;„i = 0.5 in all cases with Axm = Az,„ = A?;,„. 
For testing the convergence of the results with respect to GGEM parameters, some simulations were also 
run with 20% extra mesh points in each of the directions corresponding to aAy,„ = 0.417. We discuss next 
an important issue in multiparticle simulations, which concerns the treatment of near singular integrals - 
these integrals arise when the gaps between the particles become small. 

In a suspension of particles subjected to shear, it is not uncommon to find particle pairs separated by a 
small gap. This is true, at least occasionally, even in suspensions with a moderate volume fraction such as 
(f) = 0.15 studied here. When the gap between a particle pair becomes small, the interparticle contributions 
between such a pair from both the single as well as the double layer integrals become nearly singular; these 
nearly singular integrals require special treatment or the simulations may diverge 6^. In the literature, 
several techniques have been proposed to address this numerically difficult problem. A few among these are 
the near singularity subtraction technique 
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65| and coordinate mapping techniques 



47 1 for the evaluation 



of the near singular integrals. These techniques will, however, require very high surface mesh resolution and 
very high order quadrature techniques for the accurate evaluation of the nearly singular integrals [16 



the cost associated with these requirements can be prohibitive. An alternative approach is the use of 
a short range repulsive force, which can prevent the formation of small gaps [l^, though this is not very 



63j . The best approach appears to be an overlap correction in an 



effective in three dimensional simulations 

auxiliary step [q^. This approach is also frequently used in suspensions of rigid particles 



present work, like in past efforts 
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In the 



63j . we not only correct the overlaps, but also maintain a minimum 



gap (ft.™) between the surfaces of two particles. This approach was also employed in our recent work 3l| . 
For simulations in the current section, we set the minimum gap parameter to a small value of ft,™ = 0.05a. 

The overlap correction procedure employed in this work involves moving a pair of overlapping particles 
apart along their line of centers like a rigid particle until the minimum gap requirement is satisfied. Trans- 
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lating the capsules like a rigid particle in this auxiliary step has the benefit that the shapes of the particles 
remain unchanged, as is the orientation of the particles with respect to the flow. In general, multiple steps 
of the correction procedure is required, as the correction of overlap between one pair could result in other 
overlaps This overlap correction step involves minimal displacement of the particles, on the order 

of /i™. Given that the volume fraction studied in this work is (p = 0.15, this procedure was rarely required - 
on an average, less than 0.1 particle pairs in the A^p — 120 particle system exhibited minimum gap violations 
at any given time. We finally remark that apart from correcting the minimum gap violations in the system, 
no special treatment was accorded to the evaluation of the near singular integrals in the present effort. 

Having discussed the procedure for controlling the minimum gap in the system, we now turn to the results 
from multiparticle simulations. All these simulations were initiated by placing the particles randomly in the 
simulation box, and then they were sheared for a total non-dimensional time of t* = 7^ = 20. We show 
some representative snapshots from 120 particle simulations in Figs. (|13b .) and ([T^b ) for A = 1 and A = 5 
capsules, respectively. It is immediately obvious from these snapshots that the more viscous capsules deform 
less and also have a smaller inclination angle with the fiow direction; both of these observations arc similar 
to observations in isolated sheared capsules [sj]. The convergence of the simulation with respect to GGEM 
parameters is demonstrated in Fig. (jl3b ). where we plot the suspension apparent intrinsic viscosity [77] for 
A = 1 capsule and Np = 120 particle system for two different mesh resolutions corresponding to aA?/™ ~ 0.5 
and aAi/m ~ 0.417. The apparent intrinsic viscosity is defined as [?/] — '^^y/ (M70)i where Y,Py is the particle 
contribution to shear stress, while fi is the suspending fluid viscosity. The particle contribution to the stress 
tensor is given by [3] 

1 \ f 

= 77 X! / i^fi^j + t^i^ ~ ^){uinj + Ujn.i)]dS, (65) 

m=l 

where the sum in the right hand side is over all the particles in the system. As can be seen in Fig. (jl3b ). 
[rj] is nearly identical for simulations run with aAym = 0.5 and aAy„i = 0.417, thereby demonstrating 
the convergence of the simulation with respect to GGEM parameters. All the remaining simulations were 
performed with aAy,„ = 0.5. The effect of viscosity ratio on the apparent intrinsic viscosity is shown in 
Fig. ([TSH.). These results represent an average over the last 10 time units in the Np = 120 particle systems, 
which have a confinement ratio of 2a/ H = 0.134. For a direct comparison, the plot also shows [rj] for a dilute 
suspension of capsules obtained from single particle simulations in the same geometry. In dilute suspensions, 
a non-monotonic variation of the viscosity with A is obvious - this behavior has been demonstrated before 

n 

in the literature [^l- However, the non-monotonicity vanishes at the non-dilute volume fraction of = 0.15. 
This indicates that the contribution to the overall stress from particle-particle interactions is a monotonically 
increasing function of A, and, in non-dilute suspensions, easily compensates the non-monotonic variation 
of the isolated particle contribution. Hence caution is warranted before extrapolating trends from dilute 
systems to non-dilute systems like blood flow. 
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1.2. Computational complexity 

We devote the remainder of this section to analyzing the overall computational complexity of our algo- 
rithm. Before presenting the timing results from the detailed numerical simulations presented above, it will 
be useful to first discuss the expected computational cost associated with various steps in the algorithm, 
and consequently the overall implementation. The first step in the solution procedure, as discussed in Sec. 
(O, involves iteratively solving for the single layer density (step 1). This is followed by the computation 
of the fluid velocity u*" at the element nodes (step 2) using the single layer density q*" computed in the 
previous step. The computational cost associated with each iteration of step 1 and that of step 2 has an 
identical optimal scaling with N — N/\ x Np^ each of which is essentially determined by the computational 
cost associated with the Stokes flow solver (GGEM) described in Sec p.2p . Therefore, for a direct cost com- 
parison with the multiparticle flow problem here, we consider an auxiliary problem involving a collection of 
N point forces; both problems then have an identical computational cost scaling with N. At this stage, it 
will also be worth pointing out that the computational complexity analysis presented here closely follows 
the corresponding analysis in PME like methods 
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4l| as the underlying ideas are fairly similar. 



The overall cost associated with the GGEM Stokes flow solver is the sum of costs associated with the 
local problem and the global problem. The cost of the local solution scales as the product of the number 
of point forces N and the number of neighbors within a distance Vcut ct~^ of each of the point forces. If 
we require the computational cost of the local problem to scale as ti ^ 0{N), then the number of neighbors 
per point force must stay constant with changing system size (meaning N here). The system size N can be 
increased in two contrasting ways: (i) by increasing volume at constant density (i.e., by increasing V while 
maintaining N/V constant, V is the system volume), and (ii) by increasing density at constant volume (i.e., 
by increasing N/V while maintaining V constant). If the system size is increased at constant density, we 
require that a (or rcut) be held constant, and if the system size is increased at constant volume, we require 
that a ^ A^^/"^ (or rcut ^ N^^^"^). This scheme for choosing a ensures that the average number of near 
neighbors per point force is independent of the system size. Hence, irrespective of how N is varied, we always 
obtain ti ^ 0{N). Next, we determine the computational cost of the global solution procedure. Before we 
proceed further, it is important to realize that the error in the global solution is essentially controlled by 
the parameter aAy„|3; see Sec. p.2.2p . Therefore, for the error in the global solution to remain of the 
same order with changing system size A^, we require that aAy,„ be held constant; this implies Ay„i ^ a~^. 
Coupling this requirement with the choices of a discussed above in different scenarios, we conclude that the 
total number of mesh points involved in the calculation of the global solution n — N^NyN^ must be varied 
proportionally to A^, i.e. n ^ N. Having determined the scaling of n, we next present the expression for the 



^the error in the local solution is set by the choice of the parameter a rcut, which is unchanged with changing TV 
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computational cost of the global solution tg as follows: 

Tl n TL 

tg - jrOiN, log TV,) + j^OiN, logiV,) + — (0(iV,) + 0{Ny log TV,)), (66) 

X z y 

where the first two terms on the right hand side denote the cost associated with the FFT operations in x 
and z directions respectively, while the last term is associated with the cost of the Chebyshev-tau solver 
in the wall normal y direction. Note that the 0{Ny) cost in the expression for the Chcbyshcv-tau solver 
is associated with the quasi-tridiagonal solve, while 0{Ny log Ny) cost is associated with the computation 
of the Chcbyshcv transforms and its inverse with FFTs. Simplifying the above expression and noting that 
7T, ^ iV, we obtain the following asymptotic scaling 

n log n - iV log N. (67) 



The overall cost per iteration of step 1 or of step 2 is therefore, 



t = ti+tgr^O{N) + 0{N log N)^N log N. (68) 

The total cost of step 1 is the cost per iteration times the number of iterations required for convergence. 
Now, if the number of iterations in step 1 is independent of N, then the computational cost of the overall 
algorithm will scale as A^logiV. On the other hand, if the number of iterations in step 1 is dependent on 
the system size, say it scales as iV^, then the computational cost of the overall algorithm will scale as: 

tr^N^+nogN. (69) 

Note that the need for iterative solution of step 1 or a related second kind integral equation is not unique 
to^he present formulation, but is a general feature of any accelerated boundary integral method with A 7^ 1 

Having determined the expected scaling of our algorithm, we next report timing results from the detailed 
numerical simulations presented in Sec. (j7.ip . All runs were performed on a single core of a eight core 
machine with a 2 GHz Intel Xeon processor running Linux. We plot the time required per stage of the two 
stage midpoint time stepping algorithm in Fig. p^ . Note that the abscissa in the plot is A'plogA'j,; the 
four data points in this plot are respectively for Np = 15, Np = 30, Np = 60 and Np = 120 as discussed 
above. One can conclude from this plot that in all cases the computational cost scales approximately as 
t NplogNp (in fact, in this case, the increase in the computational cost appears to be slower than the 
expected NplogNp). Another important feature to note in the plot is the jump in computational cost as 
one moves from a matched viscosity problem to a non-matched one, which is expected as no iterations 
are required in matched viscosity problems. The number of iterations required for convergence was found 
to be independent of the system size for problems considered here, though it was found to increase with 
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Np log(Np) 

Figure 14: Wall clock time per stage of a two stage midpoint method for various viscosity ratios A and for various system sizes. 
The four data points in the plot correspond to Np = 15, 30, 60, and 120 respectively. Note that the time is plotted against 
Np log Np. 

increasing A. The simulations with A = 2 capsules required approximately 2 iterations on an average, 
while the simulations with A = 5 capsules required approximately 3 iterations. To summarize this section, 
we note that a near perfect scaling of NlogN is obtained for both matched viscosity and non-matched 
viscosity problems. A few words of caution are necessary here, though, as at higher volume fractions and/or 
at much larger system sizes, the number of iterations for convergence is expected to become system size 
dependent. It must be emphasized here that this aspect is not specific to our implementation, but is intrinsic 
to the boundary integral equation for non-matched viscosity problems. Future work on enhancements in the 
algorithm should address this by the development of efficient preconditioners. Another obvious enhancement 
in the algorithm will be its parallelization to take advantage of cheaply available multicore processors. 

8. Conclusions 

A new accelerated boundary integral method for multiphase Stokes flow in a confined geometry was 
presented. The complexity of the method scales as O(A^logA^) for the slit geometry discussed in the present 
paper. The acceleration in the method was provided by the use of General Geometry Ewald-like (GGEM) 
method for the fast computation of the velocity and stress fields driven by a set of point forces in the geometry 
of interest. Due to non-periodic nature of the domain, an alternative boundary integral formulation was 
employed, necessitated by the requirements of the acceleration technique. An efficient methodology was 
presented to compute the resulting double and single layer integrals using the GGEM technique. The 
resulting implementation was validated with several test problems. The computational complexity of the 
algorithm was verified to be 0{N log N) with timing results from several large scale multiparticle simulations. 
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Figure 15: Schematic of the problem described in [Appendix A| The figure shows the dispersed phase with viscosity fiji and a 
particle with viscosity fiB ■ The surface of the particle is denoted by S. Also shown is the undisturbed fiow denoted by u°° . 

Appendix A. Derivation of Boundary Integral equation 

Consider a two phase flow in a specified geometry as shown in Fig. ([!]). The two fluids are respectively 
denoted by A and B in the figure with viscosity and /i^. The fluid B is assumed to be enclosed by 
an impermeable interface denoted by S, which has its own characteristic properties (e.g. drops, capsules, 
vesicles, etc). Now, consider a point xq in fluid A as shown in Fig. (|15|) . Applying the reciprocal theorem to 
the disturbance velocity in the region A denoted by u^'^ = u"* — u°° , and that due to a point force located 
at Xq, we obtain 

^.f^(xo) = — ^ / /f^(x)G,,(x,Xo)d5(x) + ^ / u^^ (x)r,,,(x,Xo)7ife(x)d^(x), (Al) 
Htt^a Js Js 

where f^^ is the hydrodynamic traction at the interface on the side of fluid A associated with the velocity 

field u^^, while G is the Green's function for the specified geometry and T is the associated stress tensor. 

Next, employing the self-adjointness property of the Green's function G, we obtain the following form of 

the above equation 

uf ^i^o) = ^ I /f^(x)G,,(xo,x)d5(x) + ^ / uf^(x)r,,fc(x,xo)nfc(x)d5(x). (A.2) 

Recall that the self-adjointness of the Green's function G implies that 

Gy(x,xo) = Gji(xo,x). (A.3) 

Next, we apply the reciprocal theorem to the undisturbed flow u°° in region B and due to a point force 
located at Xg in region A. This yields, 

= / /f °°(x) G,,(xo, x) d5(x) ~^iB I <(x) T„-fc(x, Xo) nfc(x) dS{^) (A.4) 
J s Js 
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Note that in the above equation, the normal n is pointing out of region B into region A (outward normal) . 
We next note the following relation between f"^°° and f^°° 

= — . (A.5) 

Ms fJ-A 

This follows from the continuity of the undisturbed flow across the interface, and that the internal and 
external stresses are computed with viscosity fiB and respectively for the same flow field. Using this 
relationship, we write Eq. <\AA\i as 

= --^ / //°°(x)Gj,(xo,x)d5(x) + / ^/,r(x)T^yfc(x,xo)nfc(x)d5(x). (A.6) 

Adding Eqs. (|X2|) and (|X6|) . we obtain 

87r^A^^f^(xo) = - / /^(x)G'j,(xo,x)d5(x)+/^^ / (x) r,jfc(x, Xq) nfc(x) d5(x). (A.7) 

Next, we apply reciprocal theorem to the fluid velocity in region 2 (u^) and due to a point force located 
at Xq. This yields 

= / //(x) G,,(xo, x) dS{^) -fiB f (x) ryfe(x, xo) nfe(x) dS{x.). (A.8) 
Js Js 

For interfacial flows neither the interfacial velocity nor the interface tractions on either side is known, 

though we have the following boundary conditions at any point on the interface 

u-^ = = (A.9a) 

A/ = f-^ - = -f^ (A.9b) 

which essentially implies the continuity of the velocity across the interface and that the net force on an 
element due to hydrodynamic stresses is balanced by the net force due to interfacial stresses, see Sec. 
We note that the interfacial contribution f ^ is assumed to be known by the known constitutive equation of 
the interface and its known configuration. Given the unknowns and the boundary conditions, a very common 
approach is to eliminate the unknown hydrodynamic interfacial tractions f and f ^ with the known traction 
jump A/ in Eqs. (|A.7p and (jA.Sp . This leads to the widely employed second kind integral equation for 
the unknown interfacial velocity. Here, we adopt an alternative approach. Using Eqs (jA.7p and (|A.8[) we 
instead eliminate the double layer integral to obtain 

87ruf^(xo) = -— / /^(x)G,,(xo,x)d5(x) + — / /f (x) G,,(xo, x) d5(x). (A.IO) 

MA Js MB Js 

The above can be written in the following form 

u,{^o) = ^.f (xo) + ^ / - G,,(xo,x)d5(x). (A.ll) 

46 



This equation for the vefocity can be shown to be vahd inside, outside, as well as on the boundary S. The 
drawback of this equation is that both the velocity (including interfacial velocity) and the surface tractions 
are unknown. The main advantage for our purposes here is that we have switched the pole and the field 



point of the Green's function using its sclf-adjointness property (Eq. lA.Sp . For simplifying the notation, we 
next express the operand of the Green's function by q, i.e. we define q as 



Using the above definition, we write the velocity as 



uji^o) = ^if (xo) + / g.(x) G,,(xo,x) d5(x). (A.13) 
Js 

The pressure associated with the above velocity field in the region external to S can be written as 

p(xo) =p°°(xo) / 9,(x)P^(xo,x)d5(x). (A.14) 



Js 

Using Eqs. (jA.13p and (|A.14p . we can write the stress cr^(xo) as 

afk (xo ) = <yf,?° (xo )+fiA (x) Tj,k (xo , x) dS'(x) (A. 15) 

Js 

A similar expression can be written for (T^,(xo) as shown below: 

affe(xo) = a^i°°(xo) + ^ib [ q^i^) r,»fe(xo, x) d5(x) (A.16) 

Js 

We will now take the limit of the above equations as we approach the interface from either side and then 
dot it with the normal to get the tractions on cither side of the interface f'^ and f^. We can express both 
of them using the principal value of the double layer integral along with the jump condition to obtain [s^ 

//(xo) = //°°(xo) - Attha gj (xo) + fiA nfc(xo) / q^{x.) Tj,k{^o, x) d5(x), (A.17a) 

Js 



PV 



//(xo) =//°°(xo)+47rAii3<z,(xo) + /iBnfe(xo) / g,(x) T,.fc(xo, x) d5(x), (A.17b) 

Js 

where PV implies the principal value of the improper integral when the observation point lies on the domain 
of the integration. Note that the sign of the jump condition (47r(7j(xo)) depends on the direction from which 
we approach the interface relative to the outward normal defined above, i.e. whether we approach the 
interface parallel to the normal or anti to it. 

Taking the difl:erence of Eqs. (|A.17ap from (|A.17b|) and using Eqs. (|A.9bp and (|A.5|) . we have that 

-A/,(xo) = (A-l)//°°(xo)+4^A'A(A + l)g,(xo) + (A-l)/iAnfc(xo) / (?,(x) r,,;fe(xo, x) d5(x), (A.18) 

Js 

where we have now introduced the viscosity ratio A ~ ^.B/fJ-A- Rearranging the above equation, we obtain 
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a second kind integral equation for the density of the Green's function q as follows: 



where we have defined k as 



A - 1 



(A.19) 



(A.20) 



The above equation is used to solve for the unknown density q, which upon substitution in Eq. (|A.13[) gives 
the velocity at any point in the domain, including the interface. 

Appendix B. Undisturbed flow stress 

For pressure driven flows, the stress tensor is given by 



(B.l) 



where is 





f ' 


1 





e = 


1 



















(B.2) 



Note that the velocity and pressure field in pressure driven flows is given by the following expressions 

u = 4[/o-| (l - 1) , (B.3a) 

P = — jp-x. (B.3b) 

In the above equations, Uq is the centerline velocity. The surface traction f°°(x) can be obtained at a point 
X on the surface with normal vector n(x) as 



/r(x) = a-(xK(x). 
For simple shear flows, the stress tensor is given by 



(B.4) 



(B.5) 



where 7 is the shear rate. The surface traction for this case can be obtained by substituting the stress tensor 
in the above equation in Eq. (jB.4|) . 



Appendix C. Fast Spectral Stokes Flow Solver 

Here we discuss the solution procedure for the global problem in Eq. ([TU|) for a slit geometry (Fig. [T]). 
We simplify the notation in Eq. pUj) and represent it by the following set of equations: 
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- Vp(x) + ^V2u(x) = f (x), (C.la) 

V • u(x) = 0, (C.lb) 

where p is the pressure, u = (u, v, w) is the velocity, while f (x) = {fx, fy, fz) is a known function obtained 
from the known distribution of global force densities. The above set of equations are supplemented by the 
periodic boundary conditions in x and z directions, while a Dirichlct boundary condition for the velocity is 
specified in the y direction: 



u(a;, = gi(a;, z) at y = 0, (C.2a) 
u(x, z) = g2(a;, z) a,ty = H. (C.2b) 
The velocity and pressure variables are first expanded in truncated Fourier series in x and z directions 

as: 

l=-N^/2 m=-N^/2 

Similar expressions are written for f(x), w(x), p(x), /^(x), fy{x), and /^(x) by replacing uim by vim, 
Winn Pim, fxim, fyinn and fzhn respectively in the above equation. Substituting this in equation (|C.1[) and 
employing the Galcrkin method, we obtain the following owing to the orthogonality of Fourier modes: 

^2 

~ilp~ ^i{k^ - —)u = fx, (C.4a) 

- imp- ^{k^ ~ -^)w ^ fz, (C.4c) 

— + imi() = 0, (C.4d) 
dy 

where we have dropped the subscript Im in the above equation. Next, using the last equation, we 
eliminate u, v and w in the first three equations to obtain the following equation for the pressure: 

^^k^ = ^lfx + ^+^mfz. (C.5) 
The continuity equation can be replaced by the above equation for pressure along with the boundary 



condition requiring the velocity to be divergence free 



10|,l7|, 



dij 

ilu + — + imw = Oat?/ = k y = H. (C.6) 
dy 
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Alternatively, one solves the Eq. IC.5I with the following pressure boundary condition 



p = Pi at y = 0, (C.7a) 

p ^ P2 at y ^ H, (C.7b) 

though, the pressure boundary conditions above are unknown a priori, but instead they must take a value so 
that the condition in Eq. (jC.6|) is satisfied. Grouping all the equations to be solved, we have the following 

^-e-p^ ilU + §1 + *"^/-' (C-8a) 
/^"H5 fik'^u- Up = fx, (C.8b) 

fi-^ - Hk'^w ~ imp = f;,, (C.8d) 



with the following boundary conditions 



p^Pi, u = gix, V = giy, w = giz at y 0, 
p = P2, u 522;, V = 52y, w ^ g2z at y = H, 



(C.9) 



where gi = {gix, giy, giz) and g2 = {g2x,g2y,g2z) and, as with the other quantities, " denotes discrete 
Fourier transform in x and z. 



The Kleiser-Schumann influence matrix approach involves solving three set of equations as in Eq. ([C7 



with different boundary conditions as discussed shortly. In the first set, one solves the equations in lC.SI with 
the correct boundary conditions for the velocity in Eq. (|C.9p . but with homogeneous boundary conditions 
for the pressure, i.e. 

p = at y = & y = H (C.IO) 

Each of the equations for pressure and velocity components above are solved here by expanding them in 
discrete Chebyshcv polynomials and then employing the Galerkin method to obtain equations for each of 
Chebyshev modes. The appropriate boundary conditions are satisfied by employing the tau method 0,|4^ in 
which the equations for the highest two Chebyshev modes are replaced by the boundary condition equations. 
The solution for pressure is first computed, whose value is then substituted in the equations for velocity. 
Thus, in each step, one needs to solve a Helmholtz equation using Chebyshev polynomial expansion. In this 
case, the equations for the unknown Chebyshev coefficients can be reduced to a quasi-tridiagonal matrix 
equation with the last full row being full, while the rest being in the standard tridiagonal form. Also, note 
that the equations for the even and odd Chebyshev modes are decoupled and solved separately. These 
quasi-tridiagonal systems of equations can be solved in 0{Ny) time with a direct algorithm [4^. Also, the 
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use of Chcbyshcv Gauss-Lobatto quadrature points [6|, |49| for transforming a variable between physical and 
transform space ensures that FFTs can be employed for its evaluation Therefore, for each of the Fourier 
modes {l,in), the asymptotic cost of the Chebyshev-tau solution procedure scales as Ny log Ny. We denote 
this first set of solution thus obtained for Fourier mode {l,m) by {ua,i'a-,Wa,'Pa)- 

The next two sets of equations involve solving the homogeneous version of the differential equations in 



IC.81 i.e., the right hand side of the each of the equations in lC.Sl is set to zero. Moreover, the velocity boundary 
condition for these two sets of problems are also homogeneous. The only non- homogeneous equation in these 
two problems are the pressure boundary conditions. In the first of these, the pressure boundary condition 
is the following: 

p(0) = 1 & piH) = 0, (C.ll) 
while in the second the pressure boundary condition is 

p(0) = & p{H) = 1. (C.12) 

We denote these two solutions by {ub,Vb,Wb,Pb), and (uc,Vc,Wc,Pc)- It is important to note that the latter 
two set of equations are to be solved just once at the beginning of the simulation and the results are stored. 
The overall solution for the velocity and pressure is then obtained as: 



U = Ua +PlUb +P2Uc, 

V = Va + Pli'b + P2Vc, Tl\ 
W = Wa +PlWb +P2Wc, 



P=Pa+ PlPb + P2Pc 



It is easy to see that the above solution satisfies both the differential equations as well as the boundary 
conditions. The only remaining task is therefore to determine the pressure boundary conditions pi and p2. 
This is accomplished by the requiring that the velocity be divergence free at the boundary (Eq. IC.6|) . Thus, 
one obtains the following equations for the pressure boundary conditions 



CbiO) Ce(0) 

Cb(H) C,{H) 





(C.14) 



In the above equations, Cq, Cb and Cc are the continuity expression (Eq. IC.6|) evaluated at the appro- 
priate boundary point, i.e. at y — or y ^ H. For example, Ca(0) is given by 

CaiO) = UUaiO) + ^(0) + imWaiO) (C.15) 

dy 

The 2x2 coefficient matrix in Eq. (jC.14p is known as the influence matrix and the resulting matrix equation 
is trivially solved at a cost of 0(1). The solution in Eq. (|C.13[) gives the Fourier coefficients of the pressure 
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and velocity as a function of the wall normal coordinate which is then employed to obtain the pressure 
and velocity in the physical space using inverse FFTs. 



52 



References 



Allen, M. P., Tildcsley, D. J., 1989. Computer Simulation of Liquids. Oxford University Press. 

Bagchi, P., Kalluri, R. M., 2010. Rheology of a dilute suspension of liquid-filled clastic capsules. Physical Review E 81, 
056320. 

Barthes-Biesel, D., 2011. Modeling the motion of capsules in flow. Current Opinion in Colloid & Interface Science. 16, 
3-12. 

Barthes-Biesel, D., Diaz, A., Dhenin, E., 2002. Effect of constitutive laws for two-dimensional membranes on flow-induced 
capsule deformation. Journal Of Fluid Mechanics 460, 211-222. 

Blake, J. R., 1971. A note on the image system for a stokeslet in a no-slip boundary. Proceedings of the Cambridge 
Philosophical Society 70, 303-310. 

Canuto, C, Hussaini, M. Y., Quarteroni, A., Zang, T. A., 2006. Spectral Methods Fundamentals in Single Domain. 
Springer- Verlag. 

Canuto, C, Hussaini, M. Y., Quarteroni, A., Zang, T. A., 2007. Spectral Methods Evolution to Complex Geometries and 
Applications to Fluid Dynamics. Springer- Verlag. 

Charrier, J. M., Shrivastava, S., Wu, R., 1989. Free and constrained inflation of elastic membranes in relation to thermo- 
forming - non-axisymmetric problems. Journal Of Strain Analysis For Engineering Design 24, 55-74. 

Deserno, M., Holm, C, 1998. How to mesh up ewald sums. i. a theoretical and numerical comparison of various particle 
mesh routines. Journal of Chemical Physics 109, 7678-7693. 

Deville, M., Kleiser, L., Montigny-Rannou, F., 1984. Pressure and time treatment for chebyshev spectral solution of a 
stokes problem. International Journal for Numerical Methods in Fluids 4, 1149-1163. 

Doddi, S. K., Bagchi, P., 2009. Three-dimensional computational modeling of multiple deformable cells flowing in mi- 
crovessels. Physical Review E 79, 046318. 

Elman, H. C, Silvester, D. J., Wathen, A. J., 2005. Finite elements and fast iterative solvers: with applications in 
incompressible fluid dynamics. Oxford University Press. 

Fard, A. S., Hulsen, M., Meijer, H., Famili, N., Anderson, P., 2012. Adaptive non-conformal mesh reflnement and extended 
flnite element method for viscous flow inside complex moving geometries. International Journal for Numerical Methods in 
Fluids 68 (8), 1031-1052. 

Ferziger, J. H., Peric, M., 2002. Computational methods for fluid dynamics. Springer- Verlag. 

Foss, D. R., Brady, J. F., 2000. Brownian Dynamics simulation of hard-sphere colloidal dispersions. Journal of Rheology 
44, 629-651. 

Freund, J. B., 2007. Leukocyte margination in a model microvessel. Physics of Fluids 19, 023301. 
Fung, Y. C, 1996. Biomechanics: Circulation, 2nd Edition. Springer- Verlag. 

Greengard, L., Kropinski, M. C, Mayo, A., 1996. Integral equation methods for stokes flow and isotropic elasticity in the 
plane. Journal of Computational Physics 125, 403—414. 

Greengard, L., Lee, J.-Y., 2004. Accelerating the nonuniform fast fourier transform. SIAM Review 46, 443-454. 
Greengard, L., Rokhlin, V., 1987. A fast algorithm for particle simulations. Journal of Computational Physics 73, 325-348. 
Griggs, A. J., Zinchenko, A. Z., Davis, R. H., 2007. Low-Reynolds-number motion of a deformable drop between two 
parallel plane walls. International Journal of Multiphase Flow 33, 182-206. 

Hasimoto, H., 1959. On the periodic fundamental solutions of the Stokes equations and their application to viscous flow 
past a cubic array of spheres. Journal of Fluid Mechanics 5, 317—328. 

Henniger, R., Obrist, D., Kleiser, L., 2010. High-order accurate solution of the incompressible Navier-Stokcs equations on 
massively parallel computers. Journal of Computational Physics 229, 3543-3572. 

Hernandez-Ortiz, J. P., de Pablo, J. J., Graham, M. D., 2007. Fast computation of many-particle hydrodynamic and 
electrostatic interactions in a confined geometry. Physical Review Letters 98, 140602. 

Hernandez-Ortiz, J. P., Ma, H., de Pablo, J. J., Graham, M. D., 2008. Concentration distributions during flow of conflned 
flowing polymer solutions at finite concentration: slit and grooved channel. Korea- Australia Rheology Journal 20, 143-152. 
Hughes, T. J. R., 1987. The Finite element method : linear static and dynamic finite element analysis. Prentice-Hall. 
Janssen, P. J. A., Anderson, P. D., 2007. Boundary-integral method for drop deformation between parallel plates. Physics 
of Fluids 19, 043602. 

Janssen, P. J. A., Anderson, P. D., 2008. A boundary-integral model for drop deformation between two parallel plates 
with non-unit viscosity ratio drops. Journal of Computational Physics 227, 8807-8819. 

Kennedy, M. R., Pozrikidis, C, Skalak, R., 1994. Motion and deformation of liquid drops, and the rheology of dilute 
emulsions in simple shear flow. Computers & Fluids 23, 251-278. 

Kim, S., Karrila, S. J., 2005. Microhydrodynamics: Principles and Selected Applications. Dover Publicatons. 

Kumar, A., Graham, M. D., 2011. Segregation by membrane rigidity in flowing binary suspensions of elastic capsules. 

Physical Review E 84, 066316. 

Kumar, A., Higdon, J. J. L., 2010. Origins of the anomalous stress behavior in charged colloidal suspensions under shear. 
Physical Review E 82, 051401. 

Kumar, A., Higdon, J. J. L., 2011. Dynamics of the orientation behavior and its connection with rheology in sheared 
non-brownian suspensions of anisotropic dicoUoidal particles. Journal of Rheology 55, 581-626. 

Kumar, A., Higdon, J. J. L., 2011. Particle mesh ewald stokesian dynamics simulations for suspensions of non-spherical 
particles. Journal of Fluid Mechanics 675, 297-335. 

Lac, E., Barthes-Biesel, D., Pelekasis, N. A., Tsamopoulos, J., 2004. Spherical capsules in three-dimensional unbounded 
Stokes flows: effect of the membrane constitutive law and onset of buckling. Journal of Fluid Mechanics 516, 303-334. 



53 



Lac, E., Morel, A., Barthes-Bicscl, D., 2007. Hydrodynamic interaction between two identical capsules in simple shear 
flow. Journal of Fluid Mechanics 573, 149-169. 

Lambert, J. D., 1997. Numerical Methods for Ordinary Differential Systems: The Initial Value Problem. John Wiley & 
Sons Ltd. 

Le, D.-V., Tan, Z., 2010. Large deformation of liquid capsules enclosed by thin shells immersed in the fluid. Journal of 
Computational Physics 229, 4097-4116. 

Li, X., Sarkar, K., 2008. Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by 
an elastic membrane. Journal of Computational Physics 227, 4998-5018. 

Lindbo, D., Tornberg, A.-K., 2010. Spectrally accurate fast summation for periodic stokes potentials. Journal of Compu- 
tational Physics 229, 8994-9010. 

Lindbo, D., Tornberg, A.-K., 2011. Spectral accuracy in fast ewald-based methods for particle simulations. Journal of 
Computational Physics 230, 8744-8761. 

Liron, N., Mochon, S., 1976. Stokes flow for a stokcslet between two parallel fiat plates. Journal of Engineering Mathematics 
10, 287-303. 

Loewenbcrg, M., Hinch, E. J., 1996. Numerical simulation of a concentrated emulsion in shear flow. Journal of Fluid 
Mechanics 321, 395-419. 

MacMeccan, R. M., Clausen, J. R., Neitzel, G. P., Aidun, C. K., 2009. Simulating deformable particle suspensions using 
a coupled lattice-Boltzmann and finite-element method. Journal of Fluid Mechanics 618, 13-39. 

Meng, Q., Higdon, J. J. L., 2008. Large scale dynamic simulation of plate-like particle suspensions, part ii: Brownian 
simulation. Journal of Rheology 52, 37—65. 

Metsi, E., 2000. Large scale simulation of bidisperse emulsions and foams. Ph.D. thesis. University of Illinois at Urbana- 
Champaign. 

Muldowney, G., Higdon, J. J. L., 1995. A spectral boundary element approach to three-dimensional Stokes flow. Journal 
of Fluid Mechanics 298, 167-192. 

Noguchi, H., Gompper, G., 2005. Shape transitions of fluid vesicles and red blood cells in capillary flows. Proceedings of 
the National Academy of Sciences, USA 102, 14159-14164. 

Peyret, R., 2002. Spectral Methods for Incompressible Viscous Flow. Springer- Verlag. 

Pozrikidis, C, 1992. Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press. 
Pranay, P., Anekal, S. G., Hernandez-Ortiz, J. P., Graham, M. D., 2010. Pair collisions of fluid-filled clastic capsules in 
shear flow: effects of membrane properties and polymer additives. Physics of Fluids 22, 123103. 

Rahimian, A., Veerapaneni, S. K., Biros, G., 2010. Dynamic simulation of locally inextensible vesicles suspended in an 
arbitrary two-dimensional domain, a boundary integral method. Journal of Computational Physics 229, 6466—6484. 
Rallison, J. M., 1981. A numerical study of the deformation and burst of a viscous drop in general shear flows. Journal of 
Fluid Mechanics 109, 465-482. 

Ramanujan, S., Pozrikidis, C, 1998. Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: 
large deformations and the effect of fluid viscosities. Journal of Fluid Mechanics 361, 117-143. 
Saad, Y., 2003. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics. 
Silvester, D. J., Wathen, A. J., 1994. Fast iterative solution of stabilised Stokes systems part ii: Using general block 
preconditioners. SIAM Journal on Numerical Analysis 31, 1352-1367. 

Smart, J. R. , Leighton, D. T., 1991. Measurement of the drift of a droplet due to the presence of a plane. Physics of Fluids 
A 3, 21-28. 

Staben, M. E., Zinchenko, A. Z., Davis, R. H., 2003. Motion of a particle between two parallel plane walls in low-reynolds- 
number Poiseuillc flow. Physics of Fluids 15, 1711-1733. 

Stone, H. A., Stroock, A. D., Ajdari, A., 2004. Engineering flows in small devices: Microfluidics toward a lab-on-a-chip. 
Annual Review of Fluid Mechanics 36, 381—411. 

Swan, J. W., Brady, J. F., 2011. The hydrodynamics of confined dispersions. Journal of Fluid Mechanics 687, 254—299. 
Thibault, J. C, Senoeak, I., 2009. Cuda implementation of a Navier-Stokes solver on multi-gpu desktop platforms for 
incompressible flows. In: 47th AIAA Aerospace Sciences Meeting. 

Veerapaneni, S. K., Rahimian, A., Biros, G., Zorin, D., 2011. A fast algorithm for simulating vesicle flows in three 
dimensions. Journal of Computational Physics 230, 5610-5634. 

Zhao, H., Isfahani, A. H. G., Olson, L. N., Freund, J. B., 2010. A spectral boundary integral method for flowing blood 
cells. Journal of Computational Physics 229, 3726-3744. 

Zinchenko, A. Z., Davis, R. H., 2000. An efficient algorithm for hydrodynamical interaction of many deformable drops. 
Journal of Computational Physics 157, 539-587. 

Zinchenko, A. Z., Davis, R. H., 2002. Shear flow of highly concentrated emulsions of deformable drops by numerical 
simulations. Journal of Fluid Mechanics 455, 21-62. 



54 



